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Abstract. We present a numerical method to solve the equa- 
""xl" tions of general relativistic hydrodynamics in a given external 
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gravitational field. The method is based on a generalization 

of Roe's approximate Riemann solver for the non relativistic 

Euler equations in Cartesian coordinates. The new method is 

JjT applied to a set of standard test problems for general relativistic 

h-^r hydrodynamics, and is shown to perform well in comparison 

£-*> to existing numerical schemes. In contrast to existing explicit 

I/") methods the present method can cope with strong relativistic 

,-— ' shocks. By-products are: the characteristic form of the general 

relativistic Euler equations, a numerical method for special rel- 

NO ativity that can deal with strong discontinuities, a numerical 

J~^ scheme for the integration of the Euler equations in an arbitrary 

^_, coordinate system, possibly under the influence of (external) 

-_l gravity, and a novel method to incorporate source terms in nu- 

^- merical schemes. 
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0_ 1. Introduction 

03 

Astronomical jets are observed to be quite common. They con- 
sist of narrow plasma beams that transport energy from a re- 
gion near a central compact object out to distances of a hundred 
thousand or more times the dimensions of the central object. 
In active galaxies the central object is likely to be a black hole 
and energy is transported from the inner lOpc out to lOOkpc 
or more (Begelman et al. 1984). In the most extreme cases ve- 
locities are inferred that are in excess of 99% of the speed of 
light c. In galactic sources like SS433 the central object is likely 
to be a neutron star. Here energy is transported all the way to 
the extended radio source W50 (Vermeulen 1989), and bulk 
velocities of 0.26c are reached (Margon 1984). Hence there are 
both galactic and extragalactic objects in which the velocities 
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are close enough to the speed of light that relativistic effects 
should be taken into account. 

Strong gravitational fields may be important for jet forma- 
tion and propagation close to the central object. Therefore, gen- 
eral relativistic effects must be considered if the central region 
is to be studied properly. Apart from jets there are various other 
astrophysical situations, including supernova collapse, where 
it is desirable to have solutions of the equations of general 
relativistic hydrodynamics. Such solutions cannot be found an- 
alytically for most initial and boundary conditions. Therefore, 
one has to turn to a numerical approximation to the desired 
solution. Some of the cases for which analytical solutions do 
exist are used as test problems to see how accurate a given nu- 
merical method is. An overview of relativistic hydrodynamics 
was written by Taub (1978). 

From a mathematical point of view it is interesting to see 
how well the numerical concepts that were developed for the 
classical Euler equations apply to the special relativistic Euler 
equations and to the even more complex general relativistic 
Euler equations. 

Several numerical methods for relativistic hydrodynamics 
exist, e.g. the Lagrangian finite difference method of Van Riper 
(1979), and the spectral element method of Van Putten (1992). 
An Eulerian finite difference method based on upwind fluxes 
was developed and used by Hawley, Smarr & Wilson (1984a, 
1984b), henceforth HSW1 and HSW2. The main shortcoming 
of these methods is that they cannot deal with ultra relativis- 
tic shocks (HSW2). At some point it looked that no explicit 
numerical method would ever be able to deal with ultra rel- 
ativistic shocks (Norman & Winkler 1986). However, major 
progress has been made by using schemes that use the conser- 
vative form of the equations of relativistic hydrodynamics and 
use 'Riemann solvers' to follow the evolution of the flow. Ex- 
amples of this are to be found in Marti et al. (1991), Marquina 
et al. (1992), and Schneider et al. (1993). The earliest example 
was Eulderink (1988), which forms the basis of the present arti- 
cle. As we will show, this explicit method can deal with shocks 
in flows with speeds of up to 0.99999c. Also, unlike the other 
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relativistic Riemann solvers mentioned above, we consider the 
general relativistic Euler equations. 

Comparative tests in the non relativistic domain have shown 
that numerical methods that use characteristics are superior to 
other methods for violent flows (Woodward & Colella 1984). 
Experience has demonstrated that Roe's (1981b, 1986a) nu- 
merical method is particularly suited to treat strong shocks in 
the classical Euler equations (e.g. Chang & Liou 1988). Fur- 
thermore, this method stays close enough to the physics of the 
Euler equations to allow generalization to other equations that 
are physical generalizations of the classical Euler equations in 
Cartesian coordinates. The method makes no use of arbitrary 
parameters that must be tuned by experiment and that may not 
perform quite as well in the generalized situation. The arti- 
ficial viscosity introduced in some methods (e.g. HSW1 and 
HSW2) to smooth unphysical oscillations near discontinuities 
is a typical example of an arbitrary parameter that is not straight- 
forward to tune in the relativistic domain (HSW2). For these 
reasons Roe's method was taken as a starting point for the 
present method. 

An important advantage of studying the general relativistic 
equations is that their formulation in terms of a metric is so gen- 
eral that two interesting limits may be taken. By inserting the 
Minkowski metric a new method for solving the equations of 
special relativistic hydrodynamics results; by taking all speeds 
to be small compared to the speed of light one obtains a method 
for the classical equations of fluid motion in an arbitrary frame 
of reference. This is important because three-dimensional hy- 
drodynamics is very computer intensive by today's standards. 
Hence many problems can only be studied in coordinate frames 
in which the symmetry of the problem can be used to reduce 
its dimensionality. Usually this is achieved at the expense of 
a non-Cartesian coordinate frame. The present method can be 
applied in these situations. 

Because it is extremely difficult to decide how close the 
numerical solution of a complicated flow problem is to real- 
ity, numerical methods are usually tested on easier problems to 
which an exact solution is known. To demonstrate the capabili- 
ties of the method presented below, it is applied to the set of test 
problems for general relativistic codes listed in HSW1. These 
tests show that the present method overcomes said problems at 
ultra relativistic speeds. For the one-dimensional test problems 
of HSW it is superior to the method of HSW at a given reso- 
lution and time step. However HSW's method is likely to be 
faster. Comparisons based on equal CPU time have not been 
made because HSW do not mention the CPU times they re- 
quire, and because the CPU time depends on the machine, the 
implementation, and the dimensionality of the problem. Re- 
sults of astrophysical applications of our method were reported 
by Mellema et al. (1991), Icke et al. (1992) and Eulderink & 
Mellema (1994). 

This paper is organized as follows: in Sect. 2 the equations 
are introduced in the form in which they will be studied. In 
Sect. 3 some basic numerical concepts are introduced. Section 
4 describes the Roe solver for the classical Euler equation in 
Cartesian coordinates. Section 5 is about ways to determine the 



primitive variables from the conserved state vector. Section 6 
deals with the interpretation of data within a computational cell. 
Section 7 explains the analytical concept of characteristics in 
the presence of a non-constant metric. Section 8 contains the 
step by step recipe of the present method. Section 9 describes 
the function that chooses between a discontinuous and a smooth 
interpretation of the data in two neighbouring cells. In Sect. 10 
the general relativistic Roe solver is derived. Section 1 1 shows 
limiting cases of that Roe solver, such as the special relativistic 
Roe solver and the classical Roe solver in non-Cartesian coor- 
dinates. In Sect. 12 the method is compared to both analytical 
results and to other numerical methods. Section 13 contains the 
discussion and conclusions. 

In the rest of this paper the abbreviations C for Cartesian, NR 
for non relativistic, SR for special relativistic, GR for general 
relativistic and H for hydrodynamics are used. 



2. Equations of relativistic hydrodynamics 

The familiar equations of general relativistic fluid dynamics in 
the usual notation with the speed of light in vacuo equal to unity 
(c = 1), the gravitational constant times the central mass equal 
to unity (GM = 1) and the Einstein summation convention are 
(Misner et al. 1973) 



(pu a );a = 0, 

T a ?=S a , 



(2.1) 

(2.2) 



where ; denotes covariant differentiation, p is the fluid density in 
the comoving frame, u a the four- velocity of the fluid, T af3 the 
stress-energy tensor and S a the source term. Since the Einstein 
equations are divergence free, total energy and momentum of 
a system are conserved. Hence just like in the NR case the 
source term can be expressed as the negative divergence of 
a tensor. The energy and momentum supplied to the fluid by 
radiation for example, can both be expressed as the integral 
of the absorbed and scattered radiation and as the negative 
divergence of the radiation energy-momentum tensor (Mihalas 
& Mihalas 1984, p. 151 and p. 430). For use in numerical 
hydrodynamics the latter form is usually undesirable because 
an extra differentiation is introduced. 

In a space-like metric g a p normalization of the four- velocity 
can be expressed by 



3 a/3 M a M /3 



-1. 



(2.3) 



Upon writing g for the determinant of the covariant metric 
gap and using the identities (86.9) and (86.10) of Landau and 
Lifshitz (1971) one can rewrite the equations in the form that is 
used for numerical study below: 



(^g-pu a ) ia = 

and 



(V^r **),, 



s° 



(2.4) 



(2.5) 
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r^T^). 



Here T2 denotes a Christoffel symbol, defined as 

^7 — 2$ a (95/3,7 +g-y5,p + 9pi,s)- 



A more compact form of Eqs. (2.4-6) is 
in which S is the five-dimensional vector 



S = (0,S a Y ^O.S^S'.S^S*) 



3nT 



(2.6) 



(2.7) 



(2.8) 



(2.9) 



Each of the four F@ is a five-dimensional vector defined by 



pP 



V,^ , r pl ,^ 2 , T" 3 ) T .(2.10) 



The superscript T indicates transposition, to remind the reader 
that these vectors are to be used as column vectors under left 
multiplication by a matrix. With this understanding, we will 
henceforth omit this superscript. 

In this form the equations are used in the remainder of this 
paper. Here F° is called the state vector. Each five-dimensional 
vector F l is a flux in the coordinate direction x 1 . The contravari- 
ant coordinates are denoted by x a , and a; is loosely referred 
to as time. In what follows variables that are given functions 
of the coordinates (such as the metric) are collectively called 
the 'known functions' . By contrast the variables that have to be 
found by integrating the differential equation are called the 'un- 
knowns'. The source vector S contains terms due to curvature 
and, possibly, "external" source terms due to e.g. radiation. No- 
tice that there is a difference between the mathematical concept 
of a source term S and the physical concept S a : even if there is 
physically no source of energy-momentum (S a = 0) there may 
be a mathematical source term due to curvature. 

Next we make specific choices for the stress-energy tensor 
and the source term. In most astrophysical problems one is inter- 
ested in length scales that are much larger than the characteristic 
length scale of dissipative effects. Under those circumstances 
viscosity and thermal conduction may be neglected. This also 
holds if shocks are present, if one is interested in the effect of 
the shocks on the flow rather than in the precise shock structure 
(Landau & Lifshitz 1959 §87). Henceforth thermal conduction 
and viscosity are neglected in this paper. The stress-energy 
tensor is then that of an ideal fluid, i.e. a fluid with constant 
entropy along a streamline if no shocks and source terms are of 
influence: 



T af3 =phu a u 13 +pg 



a/3 



(2.11) 



where p is the pressure in the comoving frame and h the rela- 
tivistic enthalpy (i.e. including rest mass energy) in the comov- 
ing frame. The entropy of the fluid may increase due to shocks 
or due to the source term. An example is the source term that 
arises in case of the interaction of a fluid with a known radiation 
field (Mihalas & Mihalas 1984, p. 151 and p. 430). In that case 
entropy is conserved only if all radiation is scattered and none 



is absorbed by the fluid (Lindquist 1966). In general entropy 
is generated unless u a S a = 0, as is physically obvious. In this 
paper the external source term S a is introduced only for the 
sake of completeness. It is taken to be zero in what follows. 

Equations (2.8) and (2.1 1) need to be complemented by an 
equation of state. For a polytropic gas with adiabatic exponent 
T, it is given by 



ph = p + 



r 



(2.12) 



(e.g. Courant & Friedrichs 1967, Sec. 4). 

Equation (2.8) is a hyperbolic set of partial differential equa- 
tions. This is the mathematical manifestation of the fact that 
information takes time to spread in space. In this case this oc- 
curs through motion of the fluid itself and through sound waves 
(see Courant & Friedrichs 1967). An immediate consequence 
of a finite speed of communication is the possibility of discon- 
tinuities in the solution. In fluid dynamics, these often have a 
physical meaning. They are not only introduced through the 
initial state of a fluid but can also be created from smoother 
flows by nonlinearities. Of course real fluids behave smoothly. 
However the scale over which viscosity and thermal conduc- 
tivity smooth the fluid is so small compared to the scales of the 
flow in many (astrophysical) applications that a mathematical 
representation through discontinuities is accurate. 

The possible presence of discontinuities is also the reason 
why we use a form of the equations as the basis for numeri- 
cal integration which differs from HSW1: there exist infinitely 
many forms of the equations of hydrodynamics that are math- 
ematically equivalent to Eq. (2.8) . The form chosen here is 
special in the sense that it has no derivatives of dependent fluid 
variables in the source term. Thus the source term of this form 
will remain finite everywhere even if discontinuities are present 
in the solution. This is of course an advantage for numerical in- 
tegration. This form is uniquely determined in the sense that 
the combinations in which the unknowns enter are fixed. The 
state of this special form is called the 'conserved state' . As we 
shall see below, grouping the equations with respect to partial 
derivatives is essential for the central element in the present 
method: the Roe solver. 



3. Numerical approximation and operator splitting 

In this section some of the numerical concepts that are important 
for what follows are introduced. There exists a vast amount of 
literature on hydrodynamics (e.g. Courant & Friedrichs 1967) 
and on how to find a numerical approximation to solutions (e.g. 
HSW2; Roe 1986a). For completeness, and to show how the 
equations of general relativistic hydrodynamics fit into the the- 
ory, some of the numerical aspects of the problem are mentioned 
briefly. 

The most important difficulty of solving hyperbolic partial 
differential equations like Eq. (2.8) is that they allow discon- 
tinuous, "weak" solutions. In the presence of discontinuities 
the partial differential equation is ill defined, because two or 



Frits Eulderink, Garrelt Mellema: General relativistic hydrodynamics with a Roe solver 



more of the partial derivatives are infinite at a discontinuity 
in such a way that they cancel. Therefore one needs to take a 
step backwards and examine the physical situation that led to 
Eq. (2.8). 

Consider a volume Az'Aa^Aa: 3 centred around (a: 1 , x 2 , x 3 ) 
in space. The total mass, momentum and energy in this volume 
change in a time Ax° due to transport across the volume's 
boundaries and due to the source term: 



[[f (x° + Ax°,x\x 2 ,x 3 ) 

-F (x°,x l ,x 2 ,x 3 )] 
[\f\x°,x 1 + ±Ax\x 2 ,x 3 ) 

- F l (x°,x l - iAa; 1 ,^ 2 ,^ 3 )] 
f\F 2 (x ,x\x 2 + ±Ax 2 ,x 3 ) 

- F 2 (x° , x l , x 2 - ±Ax 2 ,x 3 )] 

/['V.-VV-iA.') 

- F 3 (x°,x l ,x 2 ,x 3 - jAx 3 )] 



?0/„0 „1 „2 „3J 7,1 .2 ,3, 



-l/„0 1 1a 1 „2 V j„0 j„2 . 3, 



?2,0 1 „2 1 A „2 3j t„0 , 1 ,3, 



p3/„0 1 2 3 1a 3J , . 1 7,2 



„0 j„1 j„2 j„3 



dx 3 = 0. 



This more fundamental equation may also be found by integrat- 
ing Eq. (2.8) over spacetime. In this form the equation remains 
valid even in the presence of discontinuities. This suggests the 
following numerical approximation: cover a region of space- 
time by a grid of volumes whose centres are fixed in time. For 
simplicity, the grid is taken to be equidistant in the coordinates 
([£ a ]fc = [x a ]o+kAx a ) in this paper. Then for the cell averages 
(F°) of the state one has the following second order accurate 
approximation: 

AX ° (\Fh n+i \Fh n+ " 



[(F°m\ = i(F o K Jlk 



A^ ([F2f , 



i,j+\,k 



[F \j-k* 



-i3n n " t ' 



^3l«+2 



(3.2) 



Ax 2 
Ax Q 
~ 'Ax 3 V ^hiMi ~ ^ ^,J,* 

where Q™ • k denotes the numerical approximation to the quan- 
tity Q([x°] n , [x l ]i, [x 2 ]j , [a; 3 ]j;). This formulation has an im- 
portant advantage as a numerical approximation because, just 
like Eq. (2.8), it assures that the flux that leaves one cell will 
enter another. Numerical schemes with this property are called 
'conservative'. Equation (3.2) shows also why it is advanta- 
geous to have a finite source term. The difference between the 
solutions of this numerical formulation and the true solution is 
due to the approximation of the integrals of the flux and the 
source term. 

What remains to be done to obtain a second order accurate 
numerical (integration) method is to find expressions for the 



flux and the source term which, in regions of smooth flow, dif- 
fer from the true flux and source term at the indicated points in 
spacetime by no more than terms of second order in A . Unfortu- 
nately, to determine these terms all at once is so complicated that 
almost all existing more-dimensional numerical hydrodynam- 
ics is based on evaluating these terms using one-dimensional 
simplified equations. This is perhaps the most limiting assump- 
tion for more-dimensional flow at present (Roe 1986b). How- 
ever with the first genuinely more-dimensional methods for 
CNRH just emerging (Powell & Van Leer 1989), these im- 
provements are not yet available for incorporation in a GRH 
scheme and are not included in the present paper. Instead the 
multi-dimensional partial differential equation is broken into 
one-dimensional parts. To reduce a complicated partial differ- 
ential equation like Eq. (2.8) to a series of simpler equations 
one may proceed as follows: rewrite the differential equation in 
a form where the partial derivative of the state with respect to 
time equals the sum of two operators 0\ and 2 . 



F% = (O. + 2 )(F°). 



(3.3) 



A second order accurate approximation to the solution of the 
original equation at a time step later is then found by a series 
of second order accurate integrations of the simpler equations. 
First integrate 



(3 - 1} F° = O 1 (F°) 



(3.4) 



over a half time step. Then use the resultant state as a start value 
for an integration over a full time step of 



F° = 2 (F°). 



(3.5) 



The result of this integration is then used as a start value for 
another integration over a half time step of Eq. (3.4) . If these 
equations are still too complicated for second order accurate 
integration, the above process can be repeated on the newly 
formed equations. Each next operator works on the result of the 
previous operator. In this way any complicated equation can be 
broken into ever simpler parts (Strang 1968). 

Here we apply this operator splitting to reduce Eq. (2.8) 
to a series of equations in one spatial dimension. As a first 
operator, part of the source term T is split off. Then three 
spatial operators are constructed that are each a combination of 
the flux in one direction and of a part of the source term. The 
sequence in which the operators are applied should be varied 
per time step to avoid unphysical drifts as much as possible 
(Strang 1968). Equation (2.8) can thus be reduced to one set of 
ordinary differential equations and one set of partial differential 
equations in one spatial dimension, for each spatial direction: 

F° = T, (3.6) 

F° + F\=X, (3.7) 

Y, (3.8) 

(3.9) 



F% + F 2 



F% + F\ = Z 



Unfortunately the available information is not sufficient to 
yield a unique splitting of the source term. In principle the 
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source terms may be chosen freely, subject only to the con- 
straint that T + X + Y + Z = S. This is the reason why a 
four-dimensional equation is much more difficult to solve than 
four one-dimensional equations. Each splitting of the source 
term defines a possible second order accurate update. In phys- 
ical terms, operator splitting as applied above assumes that all 
information exchange between cells takes place perpendicular 
to cell interfaces only. 

It is advantageous to make as much use of a priori infor- 
mation about the flow as possible to constrain the source term 
splitting. Which splitting is best depends on what is known 
about the flow one is trying to simulate. A form that is often 
used (e.g. HSW2) is 

T = S, 

X = Y = Z = 0. (3.10) 

A special case of such a splitting is an obvious choice for 
applications without source terms (S = 0). In that case one 
usually chooses T = X = Y = Z=Q. Some possibilities are 
discussed in Sect. 6 and Appendix C. 

Because the solution method for the three spatial operators 
is similar, the solutions of Eq. (3.8) and (3.9) are not shown 
explicitly in this paper. Their solutions follow from symmetry 
and the solution method for Eq. (3.7) . 

If no special use is made of X, for example if equation (3.7) 
is split once more to obtain an equation without source term and 
another equation without flux term, it is more cost-effective to 
choose X = and T = S. However, below the source term and 
flux in a particular direction are used in combination. Equation 
(3.7) may be rewritten in the form 



F° + (f 1 - J XdxA = 0, 



(3.11) 



where x\ is an arbitrary fixed x l value. We consider those 
special solutions of this equation which obey 



F\ = X, 

see Sect. 6 and Appendix C. 

Analogous to Eq. (3.2) , Eq. (3.11) suggests 

[(F°)ir l = [(F°w — 



(3.12) 



Ax 1 



[F ^-lAx\X^} 

[F^ + ±Ax\xf;:i\). 

(3-13) 



In this expression as in the rest of this paper, the non-varying 
subscripts (j, k) are omitted for ease of notation. The second 
order accurate approximations to the flux and the source term 
of Eq. (3.13) and to the source term of Eq. (3.6) are all that 
remain to be determined. 

It follows from Eq. (3.13) that with the possible exception 
of the very first time step, only the average state is given within 
each cell. An approximation to the flow within a cell is quite 
commonly found by assuming that the state is equal to the 



known average value everywhere in the cell or that it varies 
linearly. 

The assumption about the flow within a cell, given the aver- 
age state but not any information about the flow in neighbouring 
cells, is called the subgrid model in this paper. It is discussed 
in Sect. 6. Subgrid models do not take into account the dis- 
continuities that can be present in the global flow, but local 
discontinuities, due to state differences between adjacent cells, 
are supposed to lie on cell interfaces whenever the subgrid 
model is applied. Consequently, an approximation to the inter- 
action of flows in adjacent cells is sought that remains valid 
even if the flows in the two cells involved differ substantially. 

Mathematically, the key property of hyperbolic equations 
is that local changes to variables make themselves known to 
their surroundings at a finite speed: because of this property, 
the changes in state in cells not immediately adjacent to a cer- 
tain interface cannot influence the flux at that interface for 
some time. Thus if a small enough time step is chosen (the 
Courant-Friedrichs-Lewy [CFL] condition) the interface flux is 
independent of the initial condition outside the two cells that 
share the interface. If initially the state within each of these 
two cells is stationary, the interface flux is constant in time 
(if the metric is time independent). This flux can be computed 
from the condition that initially stationary states exist in the 
two half spaces on either side of the interface (see Courant & 
Friedrichs 1967). These states immediately left and right of the 
cell boundary under consideration are denoted by F® and F# 
respectively. The above simplified problem is called the Rie- 
mann problem. Its usefulness resides in the fact that the initial 
condition is a stationary solution away from the interface and 
that if a small enough time step is taken the information of states 
of cells farther away has not yet reached the cell boundary under 
consideration. 

The possible presence of discontinuities means that approx- 
imations to the solution of the Riemann problem based on for 
example Taylor series do not always yield accurate results. In- 
stead it is crucial that information is drawn from the proper 
places, because the state may differ drastically from place to 
place. To make this statement a bit more precise the concept of 
characteristics is extremely useful: for one-dimensional prob- 
lems characteristics are the paths in spacetime along which 
discontinuities of derivatives propagate, and which serve to 
bound the domains of dependence and influence of a given 
event, see Fig. 1. For multi-dimensional problems, characteris- 
tic surfaces exist whose properties are rather more complicated. 
Avoidance of these complications is a strong motive for the 
simplifying assumption that communication takes place nor- 
mal to cell boundaries. This allows one to work with quasi 
one-dimensional equations. 

The two-dimensional partial differential Eq. (3.11) can be 
rewritten as a series of ordinary differential equations (the char- 
acteristic equations) along the corresponding characteristic (see 
e.g. Courant & Friedrichs 1967). For ordinary differential equa- 
tions no discontinuous solutions exist. Therefore the solution of 
these characteristic equations varies smoothly along the char- 
acteristic and interpolation is useful. This property makes char- 
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Fig. 1. Characteristics in relativistic flow. At the interfaces and at the 
intermediate time-level in spacetime the three characteristic fluid ve- 
locities (forward and backward sound speed relative to the fluid veloc- 
ity and the fluid velocity itself) and the light cone are drawn. Along 
each characteristic of the linearized problem a particular combination 
of dependent variables remains constant. Hence the appropriate inter- 
face states may be reconstructed by tracing the characteristics back to 
the present time. The light cone is the boundary of influence by any 
macroscopical physical process 

acteristics extremely apt to describe flow with strong shocks. 
Schemes that exploit this notion are called characteristic-based 
schemes. 

With the help of characteristics the Riemann problem can 
be reduced to a set of nonlinear algebraical equations. How- 
ever, because the solution is needed only to find the flux in 
a situation which already was assumed to obey an idealized 
initial condition (the subgrid model), an easier "approximate 
Riemann problem" may be solved instead. Two popular meth- 
ods that solve an approximate Riemann problem are the Roe 
solver and the ENO scheme (Harten et al. 1987) based on the 
work of Osher and Salomon (1982). For strong shocks the Roe 
solver seems marginally better (Chang & Liou 1988). For that 
reason the present method is based on an extension of the Roe 
solver. 

A Roe solver makes use of approximated characteristics. 
How this works in detail is discussed in the next section where 
Roe's original solver (Roe 1981b, 1986a) for CNRH is pre- 
sented and explained. 

4. Cartesian non relativistic Roe solver 

Roe's (1981b) original CNR solver is the basis for the exten- 
sions in this paper. It is applicable to the classical Euler equation 
in Cartesian coordinates. These equations are obtained by taking 
the CNR limit of Eq. (2.8), subtracting the energy correspond- 
ing to the rest mass in the usual way. They have the form (2.8) 
with S = and 

F° = (p,pu\pu 2 ,pv?,pe c ) , 

F = (pu , pu u + p, pu u , pu u , , ph c u ) , 



F = [pu , pu u , pu u +p,pu u , ph c u J 
F = [pu , pu u , pu u , pu u + p, ph c u ) 



Here u , u 2 , w 3 are the velocities in the three spatial directions, 
while e c and h c are the classical (i.e. excluding rest mass energy) 
energy density of the fluid and the classical enthalpy 



+ (u ) + (u 



1 



+ [u 



+ [u 



T-lp' 

r p 

T-lp' 



(4.2) 



In the derivation of the equations of Sect. 3 no use was 
made of particular expressions for the state or the flux. Hence 
a numerical approximation to the update for the state may be 
found from Eq. (3.13), provided an expression for the interface 
flux at the intermediate time is found. This interface flux does 
of course depend on the particular expressions for the state and 
the flux. 

The central ideas of Roe's method derive from the Riemann 
problem for a set of linear advection equations: 



W i o + G i i = W i o+AW i i=0, 



(4.3) 



where A is a constant matrix. Each of the objects W and G is a 
five-dimensional vector analogous to the F a used above. 

Let A j; be the eigenvalues of A and e k the corresponding 
right eigenvectors. The set of equations (4.3) can be decoupled 
by projection onto these eigenvectors, i.e., by determining the 
coefficients a k such that 

W = J2 a kek, (4.4) 

or equivalently the coefficients bj. such that 

G = AW = J2 b ke k , (4.5) 

which are related by 

h=^kdk- (4.6) 

The characteristic form of the set of equations is then 

dbj = along dx l = Xjdx . (4.7) 

The eigenvectors represent the different physical waves in a 
fluid, while the coefficients are the wave strengths (Riemann 
invariants) and the eigenvalues are the characteristic velocities 
of the waves. The eigenvalue \j may thus be identified with the 
j th characteristic velocity. Equation (4.7) expresses the fact that 
the corresponding component of the flux does not change along 
a characteristic in this linear case. The flux at the boundary half 
a time step later can thus be found by tracing the characteristics 
back in time from that future point in spacetime (where the flux 
is needed) to the present time. If the time step obeys the CFL 
condition, this characteristic intersects the present time within 
one of the two neighbouring cells: 



(4.1) G J = .L ME* 1 ]*! " j^kAx , [x°] n )e k 



(4.8) 
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Two possible ways to estimate the coefficients at the current 
time are step- and linear interpolation (see Sect. 9 on which 
choice should be preferred in a particular situation): 

b k ([x% + i - ±\ k Ax°, [x°] n ) « 



£[(1 +**)&*]? + 5[(l-r*)&*]? + i, 



(4.9) 



or 



b k ([x% +i - ±\ k Ax°,[x%) n 

5[(l + "*)6jb]" + i[(l-"*)6*]" + i, 

where v k = \ k Ax i] / Ax x and a k = sign(^,t). In general one 
may write 

bk([x\ + i - ±\ k Ax°, [x°] n ) « 

i[{l + (7 j. - ((7 j. - J^)^} b k ]i 

+i[{l - c^ +(«r fc - v k )ipk}bk]2*u ( 4 - n ) 
where V' is called the flux limiter (see Sect. 9). Then if 

or equivalently 

G ?+i - G i = X^ 6 * e *^' 
this interpolation yields 

G n+ 2 _ 1 /'/^n . s~in \ 



(4.12) 



(4.14) 



The update of the state is, according to Eq. (3.13), 
Ax° 



wr 



wr 



2 Ax 1 
'Y] {[1 - cr k + (a k - v k )ip k ] b k e k }™ + i 

+ ^2 {[1 + Cfc - (c* - ^jOV'J;] b k e k Yl_i 



(4.15) 



Notice that just like the exact update this numerical update 
depends only on flux differences, and not on the fluxes them- 
selves! 

In order for this method to be of practical use one must be 
able to deal with nonlinear conservation equations. For this pur- 
pose the nonlinear equations are written in the form of Eq. (4.3), 
where A{W) is no longer a constant matrix. The method can 
then be generalized if one can approximate A by some constant 
matrix A. Roe (1981b) has formulated the properties that are 
important for A: 

1 . It constitutes a linear mapping from the vector space W to 
the vector space G. 

2. As W L -H- W R -h- W : A(W L ,W R ) -► A(W), where 
A = dG/dW. 

3. For any W L , W R : MW R , W L )(W R - W L ) = G R - G L . 

4. The eigenvectors of A are linearly independent. 

Condition (1) is necessary to simplify the original equations to 
a linear set. Condition (2) is necessary to recover smoothly the 



linearized algorithm from the nonlinear one. Conditions (3) and 
(4) were proven to be necessary and sufficient for a conserva- 
tive scheme with correct jump conditions across a shock (Roe 
1981a). 

In practice property (4) can be checked afterwards. Matrices 
A with properties (1), (2), and (3) can be constructed in several 
ways (Roe 198 lb). The method that was used by Roe for CNRH 
and that is used in this paper for GRH is based on the exact form 



(4-10) A(pq)=(p)Aq + (q)Ap, 



(4.16) 



where p and q are any quantities, Ap is the difference in quantity 
p, and ( } denotes the arithmetic mean. 

Once such a matrix A is constructed for the particular set 
of equations, its eigenvectors and eigenvalues may be deter- 
mined. The projection coefficients are obtained by projecting 
the state difference (4.12) or equivalently the flux (4.13) on the 
eigenvectors. The flux and the state update are then given by 
Eqs. (4.14) and (4.15) respectively. 

The Rankine-Hugoniot jump condition for exact solutions 
states that 



(4.13) G R -G L = X(W R -W L ), 



(4.17) 



where A is the propagation velocity of a stable discontinuity. 
The key point of Roe solvers is that properties (3) and (4) guar- 
antee that the Riemann problem is linearized in such a way that 
this correct propagation velocity is found if the initial disconti- 
nuity is a pure contact discontinuity or a pure shock, even for 
nonlinear problems. By comparing Eq. (4.17) to Eqs. (4.13) and 
(4.12) it may be deduced that \j = A if only one a,j ^ 0. Thus 
for a discontinuity of only one Riemann invariant Roe' s method 
finds the correct characteristic speed. This property is also the 
reason why the Roe solver needs the conservative form of the 
equations, i.e. the form of Eq. (2.8): the exact jump conditions 
only hold for this form. 

By expressing the state and the flux in terms of a parameter 
vector Kc(l , h c ,u l , u 2 , w 3 ), where Kc = yfp, and by applying 
Eq. (4.16) , Roe (1981b) found the following explicit expres- 
sions for CNRH. The eigenvectors are 



ei 



ei = 



e3 = 



e 4 = 



es = 



where 



l,h c - s c v c ,v c - s c ,v c ,v c ), 



„2 „,3 



l,h c + s c v c ,v c +s c ,v c ,v c ,, 



\^ 2 \[vi) 2 + {vlf + {vlf\ y c y c y c >, 



1 LA) 



,1 „,2 „,3 



0,^,0, l.Oj, 

0,^,0,0, A 



(4.18) 
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s 2 c = (T-l){h c -i[(vl) 2 + (vl) 2 + (vl) 2 }}, 

, (Knvf) 



C (Kc) ' 
~ {K c h c ) 
C= (Kc) ■ 
The characteristic velocities are 
\i =vl~ s c , 
A 2 = v\ + s c , 
A 3 = A 4 = A 5 = v\. 



(4.19) 



(4.20) 



The projection coefficients of a state difference F R — F^ 

(6, A e , A 1 , A 2 , A 3 ) are 

a\ = 



a 2 



^M^ + W) 2 ^) 2 ]* 

-KA' + ^A 2 + ^A 3 -A e ]} 

^{i[K) 2 + (, 2 ) 2 + (, 3 ) 2 ]^ 



2s c 



[^A'+^ 2 A 2 + ^A 3 -A e | !> + ■ 



']}■ 



A 1 -vlS 



2s, 



F ' M h c -\(v\) 2 + (vl) 2 + (vl) 2 ]\6 



a 3 - —j— y y, c - \yv cJ , yu c 



M) 2 ]}. 



+ [^A 1 +v 2 A 2 + v 3 A 3 - A e 



a 4 = A 2 - v 2 5, 



a 5 



A 3 - u 3 6. 



(4.21) 



With Eqs. (4.2 1 ) , (4.20) and (4.18) one can derive numerical 
approximations to the interface flux (and to the update of the 
state) from Eqs. (4.14) and (4.15) if F° is substituted for the 
state W and F ' for the flux G. The flux difference at an interface 
is found and decomposed under the assumption of a constant 
state within a cell. Next a correction based on the interface 
flux difference and the upwind interface flux difference for the 
relevant characteristic is applied to the flux, by means of the 
flux limiter tp. This correction models the non-constant flux 
variation. For nonlinear equations this is not the same as first 
choosing a suitable flux model and then finding interface states 
as initial conditions for a Riemann problem, as suggested by 
e.g. Van Leer (1973, 1974, 1977a, 1977b, 1979). 

Unfortunately, a Roe solver does not always yield a physi- 
cally allowed flux. An example is encountered in the Riemann 
problem at an interface of two cells filled with a gas with adia- 
batic exponent T = 1 .5 and having pl = 1 , «l = —2, pi, = 4/3 
and pr = 4, ur = 1 , pr = 13/3 respectively. After calculation 
the first order upwind flux with the Roe solver it is found that 
there is no mass flux and yet there is an energy flux. Clearly no 
physical set of variables exists that can reproduce this Roe flux. 
An explanation is the fact that the Roe solver represents expan- 
sion waves by expansion shocks with one velocity (see Einfeldt 



et al. 1991). One might consider the use of an exact Riemann 
solver in cases where the Roe solver cannot cope but including 
special cases would make the numerical code rather slow. Yet 
for most problems the Roe solver is remarkably trouble free. 
Thus a powerful scheme for the classical Euler equations in 
Cartesian coordinates results. 

The above ideas depend on characteristics rather than on 
the precise functional form of the state and the flux. Hence Roe 
solvers may be extended to hyperbolic equations other than the 
CNR Euler equation. Before Eqs. (4.14) and (4.15) can be used 
to obtain a scheme for the GR Euler equations, the eigenvectors, 
eigenvalues and projection coefficients for the GR equations 
have to be found. This is done in Sect. 10. However, first two 
difficulties need to be addressed. The first of these also appears 
in SRH and has to do with the relativistic form of the state. 
It is discussed in Sect. 5. The other difficulty also appears in 
NRH in non-Cartesian coordinates and is caused by the use of 
curvilinear coordinates. It is discussed in Sect. 6. 



5. Determination of the primitive variables 

To find an update for the state it is necessary to determine the 
state or flux differences, the source term and the characteristics 
of the flow. All of these can be expressed explicitly in terms of 
the primitive variables p, p, m'/ u °i u 2 /u { \ m 3 /m°. However, 
only the (average of the) conserved state is known at each time 
step, as is evident from Eq. (3 . 1 3). In RH, in contrast to NRH, the 
primitive variables are not easily determined from the conserved 
quantities, because the enthalpy enters the spatial components 
of the conserved vector in RH. HSW do not encounter this 
problem because they use another (non conservative) form of 
the state. The reasons for choosing the conservative form were 
explained in Sections 2 and 4. 

Appendix A is dedicated to a discussion of how to determine 
the primitive variables given the conserved quantity vector. 
Three methods are described there: an analytical method, a one- 
dimensional Newton-Raphson method and a six-dimensional 
Newton-Kantorovich method. The one-dimensional Newton 
method was applied to obtain the results in this paper. 



6. Subgrid model 

To predict the flow on a grid of several computational cells one 
must be able to predict the interaction of the flows in neighbour- 
ing cells. This interaction depends on the flow within a cell. In 
the CNRH scheme discussed in Sect. 4 for example, two inter- 
face fluxes each based on information of just one of the two 
bounding cells are combined by a Roe solver via Eq. (4.14) to 
construct the desired flux for integration in Eq. (3.13). There 
are three equations left to solve: (3.6), (3.11) and (3.12). Of 
these Eq. (3.6) as a whole and Eq. (3.12) (the integration of 
the source terms) pertain to one cell. Only after these integra- 
tions have produced the current interface state on both sides 
of the interface can the interaction between bounding cells be 
computed. Hence, the flow is computationally approximated in 



Frits Eulderink, Garrelt Mellema: General relativistic hydrodynamics with a Roe solver 



two distinct steps; the approximation of the flow within a cell 
and the approximation of the interaction between these flows in 
neighbouring cells. This section is dedicated to the former step. 
The subgrid model describes the assumed behaviour of 
the primitive variables within a cell as a function of time for 
Eq. (3.6) and as a function of place for Eq. (3.12). 



6.1. General considerations 

The numerical methods presented in this paper aim to minimize 
complications due to the metric and source terms in calculating 
the interactions between adjacent cells by solving physically 
allowed (approximate) Riemann problems at specific points in 
spacetime. No complications due to non-physical states or dif- 
ferences in the metric enter the Roe solver. To have physically 
allowed states is important because an approximate Riemann 
solver may not be expected to yield a physical interface flux 
at the intermediate time, based on a non-physical current inter- 
face state. Any correction that has to be applied to make the 
final update physical again is arbitrary by its nature. HS W force 
their results to be physical by adapting u° ("velocity renormal- 
ization", HSW2: Eq. (11). 

The differences in the metric and the source terms are taken 
into account by subgrid models. These models are used to de- 
termine the interface states from the known average state. In 
principle they can also be used to extrapolate the calculated 
interface fluxes to the (centres of the) respective cells. Hence 
the complications due to the metric and source terms are con- 
centrated in the subgrid model and the Roe solver just needs 
to deal with localized GRH equations. It is anticipated that as 
a result of this approach the remaining ingredients of the im- 
plementation, such as the flux limiter, can stay as they are for 
CNRH. 

In CH without source terms the subgrid model is so simple 
that it hardly received any attention in Sect. 4. All variables 
were assumed to be constant within a cell (prior to decomposi- 
tion onto the eigenvectors). This model is not adequate in more 
general geometries and/or coordinates. This is especially ap- 
parent in curved spacetime: a vector at one point in spacetime 
is not necessarily a vector at another point, and a linear combi- 
nation of two vectors that exist at different points of spacetime 
generally does not yield a vector in a location in between. 

Desirable properties of subgrid models are easily identified 
by considering CNRH without source terms. There the usual 
choice is a subgrid model with constant variables given by the 
cell-averaged state. This model has four advantages: 

1 . The subgrid model is computationally cheap. The only work 
involved in this subgrid model is the determination of the 
primitive variables that belong to the given averaged state. 

2. The average of the state over the cell is equal to the pre- 
scribed cell-average. 

3. The prescribed state is a stationary solution. This ensures 
firstly that the interface problem is a Riemann problem 
and secondly that a one-dimensional stationary solution is 
preserved. 



4. The prescribed interface states are always physically al- 
lowed. 

In CSRH without source terms the constant variable model 
retains these advantages and consequently it is the obvious 
choice. However in case a non-zero source term or non-constant 
known functions (metric) enter the problem, the four properties 
are in conflict. For example, a constant state cheaply assures 
the correct cell-averaged state but is in general not a stationary 
solution. Hence, some of the desirable properties above have to 
be sacrificed. 

On the one hand one can decide to discard property (3). 
In that case the constant state model can be accepted. Even 
then extra calculations are required if the metric varies because 
the same state will produce different primitive variables and 
consequently different fluxes at different places within a cell. 
In addition, this subgrid model cannot guarantee a physical 
interface state, i.e. property (4) is compromised as well. Hence, 
the disadvantages of this subgrid model are that it diffuses even 
a one-dimensional stationary solution and that non-physical 
interface states cannot be ruled out. 



Fig. 2. Subgrid model based on extrapolation. The cell average data is 
assumed to yield the variables at the cell centre. To obtain the vari- 
ables at the interfaces the variables at the cell centre are extrapolated 
via an assumption, for example that of one-dimensional stationary 
flow. These interface variables are subsequently used as input for the 
(approximate) Riemann problems 

On the other hand, to ensure minimal damage to stationary 
one-dimensional solutions, one may use these as a basis for a 
subgrid model. These solutions are in general nonlinear func- 
tions of the state. Therefore it is no longer feasible to enforce the 
cell-averaged state to be correct. Instead the cell-averaged state 
is assumed to be representative of the state at the cell centre. 
The subgrid model is then fixed by a value rather than by an 
integral equation as condition to select the proper solution from 
those permitted by the differential equation. The subgrid model 
is now effectively an extrapolation of the current data at the cell 
centre in space or in time. Schematically this extrapolation is 
indicated in Fig. 2. 

Notice that, even though the cell-averaged state may not be 
correct, such a model yields a conservative method as long as the 
flux added to a cell is subtracted from a neighbour. A possible, 
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relatively cheap implementation of this kind of subgrid model 
is to use the primitive variables that belong to the average state 
to define a constant average source term for the flux, i.e. 



IF 1 ! 



[F% 



+ \^X- 



(6.1) 



This subgrid model preserves one-dimensional stationary so- 
lutions to second order, but only has a first order correct cell- 
averaged state, and does not ensure a physical interface state. 

Hence, both subgrid models above cannot guarantee a phys- 
ically allowed interface state. In GR, parallel transport is the 
physical way to transport vectors from one point to another. In- 
deed, parallel transport is a possible subgrid model that yields 
physical interface states. However parallel transport makes no 
use of the symmetries of the flow that usually underlie the choice 
for non-Cartesian grids. This can easily be seen for curvilinear 
coordinates in flat spacetime: in that case parallel transport 
leaves vectors constant along straight lines, independent of the 
symmetry of the coordinates. As a result a radial velocity vector 
is no longer radial after parallel transport in any other direction 
than the radial direction, even if spherical coordinates are used. 

The transport presented below yields physically allowed 
interface states and makes use of symmetries of the flow with 
respect to the chosen coordinates. The freedom inherent in the 
source term splitting is used to choose a special case that is based 
on one-dimensional stationary solutions. As a result the one- 
dimensional Eqs. (3.6) and (3.12) can be integrated analytically. 

6.2. Algebraic extrapolation 

Both the extrapolation in space and in time may be done alge- 
braically if (quasi) one-dimensional flow in curved spacetime 
is assumed, and the source terms are split according to their 
metric derivative: 

2 9"'9pj, 

Y = ^ (0, -p^gp 7 , 2 T^ + jg--g Py ,: 



T = 



Z = 



1(0,- 
1{0,- 
1{0,- 
1{0,- 



-g a/3 gp 7 fi--t 






,a/3 



gp y , 3 T* + y<"g Py , 3 Tf>i) 



(6.2) 



In the spatial directions this extrapolation will henceforth be 
called stationary extrapolation, because it is the correct subgrid 
model for one-dimensional stationary flow (i.e. the state and 
the metric are time independent). Applied in the time direction 
it is called 'homogeneous extrapolation', because it is correct 
for one-dimensional homogeneous flows (i.e. flows in which 
the flux and the metric are independent of the spatial coordi- 
nates). This splitting has the additional advantage that T = 
if go = and X = if g t \ =0. Thus if the metric is time 
independent the time extrapolation can be skipped, and if the 
metric is homogeneous the space extrapolation can be skipped. 
More important than this last property is the fact that if an- 
alytical extrapolation is used the metric is taken into account 
exactly. On a Cartesian grid an unphysical prediction occa- 
sionally happens due to approximations in the Riemann solver 
(see Sect. 4) or due to the flux limiter. In GRH however, not 



taking the metric exactly into account is an additional very im- 
portant source of these errors. Near a black hole for instance, 
some components of the metric tend to infinity. Consequently 
an arbitrarily small approximation error in the source term can 
dominate the flow close to the horizon, often resulting in the 
erroneous prediction of states that cannot physically exist. 

We now derive the equations for integration in the 1- 
direction. Extrapolation in the other directions is similar. By 
assuming that while integrating in one direction (here the 1- 
direction) all functions are constant in all other directions, it 
follows from Eq. (3.11) that for a stationary solution the mass 
flux is equal to a constant D: 



-gpu 1 = D. 



(6.3) 



Replacing the Christoffel symbols by their definition Eq. (2.7) 
and taking only the 1 -derivatives into account for the extrapo- 
lation in the 1 -direction yields 



gea,iV~gT a + g e a(V-g~T a ) i i 




9ta,W—g-L - g ea g g-ys.iV-g 1 




Jffeaff g-yS.lV-g 1 




gea,iV=gT la - fler.lV^T 1 ? =0, e 


= 0,2,3. 




(6.4) 



Hence ^/—gT\ is constant. In combination with Eq. (6.3) this 
implies for D ^ that 



(hu e ) A =0, 6 = 0,2,3, 

If D = this relation is adopted. For e = 1 one has 



(6.5) 



(V=gT\) tl = \ 9ySA ^^T^ = 



5 - —i^ s 



-gTys 



i 

2 A 

id ,n 



-9t [(g^^Ujhus)^ - g 1 \hu 1 hu s ) i i] 



-gpg-ys 

= V-gph l i + \/^gpu\hui) i i + ( v ^) i ip, (6.6) 
where in the last step the identity 



1 7<5 

29yS9 ,i = 



1 



=(V=F), 



(6.7) 



is used. Because the mass flux D is constant, one may also write 
(V=gT\) A = (£>/mi),i + (V=Fp),i 

= V-gpu l (hui) t i +(v / -3P),i- (6-8) 

Equating the two expressions (6.6) and (6.8) for ^—gT\, one 
finds 

(lnp-r 111^,1=0, (6.9) 

and hence 

p = K P r , (6.10) 

where k is a constant. This relation is the requirement that en- 
tropy should remain constant in stationary interpolation (shocks 
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are supposed to lie on cell boundaries only). This means that 
this subgrid model prescribes isentropic flow. This is not the 
case for an arbitrary splitting of the source terms. The entropy 
is constant for the splitting used here because it systematically 
ignores all changes in other directions than the direction of 
extrapolation. 

The five relations (6.3) , (6.5) and (6.10) yield a set of 
algebraic equations for the stationary solution of the (quasi) 
one-dimensional flow. To find the primitive variables at an ar- 
bitrary location, one can combine these relations into an implicit 
equation for p: 



q(p) 



hi 



h 2 



\D\ 



Q- 



9V9 

Here ho is the value of the enthalpy if u l 
D = 0) 

(g u hu e ) 2 



hi 



g ev hu e hu v , 



(6.11) 
(or equivalently 

(6.12) 



where the summation only involves e, 77 = 0, 2, 3. The metric 
is evaluated at the point where the primitive variables are to be 
found. Physical solutions should have 1 < h < ho. Given the 
density, the other primitive variables follow simply from the 
constants of the flow as 



p= up 



D 



-gp 

hu f 



1 + fW" 1 



(6.13) 
(6.14) 
(6.15) 



6.3. Solution of the algebraic equations 

The constants k, D, huo, hu2, and huj, are determined by as- 
suming that the cell averages are representative for the cell cen- 
tre. Next, Eq. (6.1 1) can be solved for p, by a Newton method 
for instance. The solution is the intersection of the left-hand side 
function q, that is zero at p = and at po > and that has a neg- 
ative second derivative in that interval, and the constant on the 
right-hand side. Hence there are either zero, one, or two roots. 
It is shown below that when two roots exist one corresponds to 
a supersonic and the other to a subsonic velocity. In case one 
root exists the flow is sonic. This situation is similar to the well 
studied problem of classical steady flow along a streamline: for 
a given mass flux, energy flux and entropy there are either zero, 
one, or two possible states. In case there are two possible states 
one corresponds to a subsonic and one to a supersonic solution. 
The necessary and sufficient condition for the existence of 
roots may be found by examining q. The maximal mass flux 
that can be reached for given hu e , e = 0, 2, 3 and k is found by 
considering the p derivative of q: 

dq P(h) 



h 2 



hi 



(6.16) 



h 2 



where 

P(h) = -h 3 + (2 - T)h 2 Q h + h 2 Q (T - 1), 
and we used the relation 



dh 



(r-i)(h-i). 



(6.17) 



(6.18) 



The polynomial P(h) has at least one root /i* in the range 
1 < h* < ho, because P(l) = h\ - 1 > and P(h ) = 
— (T — 1)/iq(/io — 1) < 0. In fact there is only one root in 
the physically allowed region because the h derivative of P 
changes sign at h± = ±y(2 — T)/3ho- This means that P has 
one extremum for positive h. Thus precisely one root /i* exists 
that has 1 < h* < ho (other roots have h < 0, if they exist). 
The root formula for cubic equations shows there are one, two 
or three roots in total, respectively if ho/ hi < 1, ho/ hi = 1 or 
ho/ hi > 1, where 



hi = 



3(r - 1) 

2(2 - r) v 2 - r 



(6.19) 



It is computationally faster not to use the cubic root formula 
to find the physical root but rather a Newton method. A good 
starting value for the Newton method is h = ho- For that value 
PP" > and both P' and P" do not change sign in the interval 
(/i* , ho), so the root is approached from one side (here h > h*), 
see Abramowitz and Stegun (1970). Hence one always finds 
the required largest root /i* . Given this /i* the corresponding 
density may be found from 



P* = 



1 



Tk 



-(K - 1) 



(6.20) 



Then the maximal mass flux that can be supported by hu e , e •■ 
0, 2, 3 and k is 



\DJ = 




(6.21) 



If the required mass flux \D\ > \D*\ then no solution exists; 
if \D\ = |£)*| one solution exists that has precisely the sonic 
velocity (see below), and if \D\ < \D* \ both a supersonic and a 
subsonic solution exist. These solutions are found by applying 
a Newton method to Eq. (6.11) with start value p = and 
P = Pmax respectively, where p max is given by 



1 + 



-KO 1 '- 1 
^Pmax 



(6.22) 



1 + 



By examining the relevant derivatives it may verified that both 
solutions are approached from one side with these initial values. 
The simpler starting value given by h = ho cannot be used 
because dq/dp is infinite for p = po, see Eq. (6.16) . This would 
cause the Newton method not to update the approximation to p, 
even though Eq. (6.1 1) is not obeyed. The starting value p m ax 
is good because it is smaller than po and greater than or equal 



12 



Frits Eulderink, Garrelt Mellema: General relativistic hydrodynamics with a Roe solver 



to the largest root of Eq. (6.1 1), it has q(p m ax) < Q and it has a 
well defined dq/dp for Q ^ 0, while for Q = it is the correct 
solution. 

When no solution exists, one has to conclude that the cho- 
sen subgrid model is false. In GRH there can be two physical 
reasons for the absence of solutions. In flow around a black 
hole, for example, the fluid may not have sufficient energy to 
reach a certain larger distance from the hole. This is expressed 
by the failure of the interpolation procedure because ho < 1 . 
To force a solution when ho < 1, ho is multiplied by the ratio 
(l+e)//io to yield /io = 1+e, withe a small positive number. This 
will yield very small interface fluxes. To avoid inconsistencies 
between the interface quantities and the central quantities, the 
central state is also replaced by a state with the original D, k 
and w e , e = 0, 2, 3 but with an enthalpy that is multiplied by the 
same ratio. Small changes to the state can be explained away 
because it was assumed that the cell-average of the state is the 
actual state at the cell centre. 

The other possible cause is that the cell's surface area is 
too small for a fluid of a fixed energy to transport the required 
mass flux. In that case ho > 1, but the mass flux D is too large 
to allow a solution. The easiest way to find an extrapolation 
anyway is to put in a sonic solution. Numerical experimenting 
showed that this can best be done by adapting the mass flow 
rate to \D\ := |-D*|. Again to avoid inconsistencies between 
the extrapolated and the original flows, the original state is 
replaced by a state with the original k and hu e , e = 0,2,3 
but with |D| = \D*\. The subsonic branch is chosen for the 
corresponding state if it is upwind of the interface and the 
supersonic branch if it is downwind. Thus the existence issue 
is settled. 

Two solutions exist in nearly all cases encountered in prac- 
tice. Equation (6.11) may be rewritten in terms of velocity as 



(6.23) 



ho J u l u l + g n 



Since the speed of sound a is given by 

2 _ i> (r - \){h - i) 



ph h 








the Mach number of the flow is 






r. 1 u l u l h 
M — = 

a 2 (u l u l + g 11 ) (T-l)(h- 


-1) 


1 - 


Uoj . 



(6.24) 



T- -(6.25) 
h oJ 

By comparing this expression to the polynomial P above it may 
be seen that M 2 > 1 if P > 0, i.e. if p < p„, and M 2 < 1 if 
P < 0, i.e. ifp > p*. Hence if two solutions of Eq. (6.17) exist, 
one is subsonic and the other supersonic. The selection of the 
proper root is the remaining issue. 

The choice between the sub- and supersonic solution re- 
quires attention. It does not suffice to adhere to the same branch 
as the central state in all cases. Under the influence of a source 
term, smooth flow can go from a subsonic to a supersonic regime 
through a sonic point. A well-known example of such a smooth 
transition occurs in the de Laval nozzle. Therefore, one must 



check whether such a sonic point is present in the flow within a 
cell. This is done by solving D* (x ) = D. After a starting value 
x 1 is chosen, D* is determined by solving Eq. (6.17), (6.20) 
and (6.21) for that a; 1 . Next x l is updated by a Newton method: 



x l :=x l 



where 



D*-D 



D 



*,i 



D 



*,i 



with 




(6.26) 



(hi) , , (6.27) 



,lef,„, \/JV 



_ 2(g^hu e )(g^hu v ) (g^hu e ) 2 u ev 



(6.28) 



If a solution x 1 is found which lies between the cell centre 
and its interface, the solution is assumed to switch branches. 
This procedure is not very sensitive to numerical accuracy, 
because flow in the neighbourhood of a point with flow on one 
branch is on the same branch unless the flow is almost sonic. 
In that case the difference between the two branches is small. 
In fact, accidentally picking the wrong branch every once in a 
while may be thought of as a stability check on the flow: if it 
makes a considerable difference, the flow pattern one is trying 
to approximate probably is not stable. 

6.4. Disadvantages of stationary extrapolation 

Stationary extrapolation also has two disadvantages. The first 
is of a general nature and the second is only a disadvantage if 
extrapolation of a flux is required for which the corresponding 
primitive variables are not known. Such a situation occurs if an 
interface Roe flux is to be extrapolated to the cell centre. 

The first disadvantage is related to a typical superresolution 
problem. An assumption is used to get more resolution than 
the grid actually contains. In case that assumption is incorrect 
a penalty is paid. A specific example is given in Appendix B. 

The second disadvantage is that for algebraic extrapolation 
the primitive variables must be known. The Roe solver only 
yields the upwind flux and not the upwind state (see Sect. 4) 
at the interface of two cells. The primitive variables can be 
determined uniquely from the state, but in general not from the 
flux. This problem is discussed in more detail in Appendix C. 

To avoid having to determine the primitive variables from 
the flux we use the same source term that was used to get 
information from the cell centre to the interface to transport the 
flux in the opposite direction. If 



[X]- 



1 



Ax 1 



[F l Cl 



[FTA 



then 



,1,-^r^lW i Ax i [X f: 



[*x * = [fx;, 



(6.29) 



(6.30) 
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and 



[F ^ = [F i ]:r + l Ax i [x t 



(6.31) 



As is discussed in appendix C one must add a small correction 
term to this source term to make the resulting integration of the 
partial differential equation second order accurate. 

6.5. Concluding remarks 

So far only extrapolation in the spatial direction has been dis- 
cussed. Most of what has been said is, mutatis mutandis, ap- 
plicable also to time extrapolation. The problems due to non- 
uniqueness do not occur, because the primitive variables are 
uniquely determined by the state. For time independent metrics 
it is trivial. A schematic overview of all extrapolations involved 
is shown in Fig. 3. 



Fig. 3. Schematic overview of the extrapolations involved to obtain an 
update. The arrows indicate extrapolation of information in one point 
of spacetime to another. The time levels n + |— and n + |+ indicate 
are just before and just after the approximate Riemann problem at time 
n + ^respectively. For time independent metrics the extrapolations in 
the time direction are trivial 

In this section a physical interface state was found by 
stationary/homogeneous extrapolation of the cell centre data. 
Hence if an exact GR Riemann solver is applied to the interface 
states, a physical interface flux is guaranteed. If this interface 
flux is extrapolated back to the cell centre (the primitive vari- 
ables are known after an exact Riemann solver) a physical 
update can be guaranteed in GR. 

If a Roe solver is applied, however, a physical update cannot 
be guaranteed because of two reasons: first, a Roe interface flux 
is itself sometimes unphysical. This problem is not specific to 
the extensions discussed in this paper. An example in CNRH 
was given in Sect. 4. Second, the Roe interface fluxes are not 
extrapolated back to the cell centre by stationary extrapolation. 
Instead we used Eqs. (6.31) and (6.30) . 

Numerical experience shows that with these subgrid models 
an unphysical update rarely occurs. If it does occur, a new 
update is calculated by using the flux that corresponds to the cell 
centre and the flux that is found after extrapolating the upwind 
state all the way to the cell centre. If this does not improve 



matters this flux difference is applied for a progressively smaller 
time step, until a physical update is found. By ensuring that the 
time step is zero after a finite number of trials a physical update 
is always found, because if all else fails the new state is kept 
equal to the present state. Admittedly, this is rather ad hoc. 

7. Approximate Riemann problem 

In this section the interaction of the fluid in adjacent cells is 
described. As stated is Sect. 3 discontinuities can be present in 
the flow. The subgrid models presented in Sect. 6 do not take 
these into account. Hence discontinuities are assumed to lie on 
cell interfaces. As a result an approximation to the interaction 
of flows in adjacent cells is sought that remains valid even if 
the flows in the two cells involved differ substantially. 

The equations of GRH are more complicated than the 
CNRH equations discussed in Sect. 4 in the sense that gen- 
erally a source term is present and that the state vector and the 
flux depend on a known continuous metric g a p . Denote these 
known continuous functions by g. Then one has 



W(u,g) i0 + G(u,g) A = B(u,g). 



(7.1) 



It suffices to find the evolution of the unknowns u since the flux 
at the cell interface at the intermediate time can be built from 
combinations of these unknowns and the known quantities. The 
initial values of the unknowns follow from the initial conserved 
quantities. It is possible to isolate the differentiated unknowns 
by rewriting Eq. (7.1) as 



dW dG 

-z— u ,o+^— u-i = J, 
du du 



where 
J = B 



dW dG 

'9,Q ~ ~H~9,i- 



(7.2) 



(7.3) 



og dg 

Let ej; be the right eigenvectors and \j. be the eigenvalues of 

dG fdW\~' 



a \ a / ' (7 - 4) 

au \ au J 

and let 

dG 

— — du—Jdx =Sej;d6j;, (7.5) 

au 

then the characteristic form of this equation is (cf. Eq.[4.7]) 

dbj = along dx l = \jdx°. (7.6) 

Notice that if the known functions are time independent, one 
has 



dG 
du 



du — Jdx = dG — Bdx , 



(7.7) 



and Eqs. (7.7) and (7.6) indicate that the interface states de- 
termined by stationary extrapolation may be used to determine 
Roe fluxes at these interfaces. A possible time dependence of 
the metric is taken into account by Eq. (3.6). 
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8. Implementation for general relativistic flow 

In Sections 3, 5, 6, and 7 the ingredients of numerical methods 
for GRH were discussed. In this section these elements are 
integrated. The update of the fluid state after a time step At is 
determined as follows: 

1. Given the state [F°] i 2 determine the primitive variables 



2. Calculate the spatial source term X { 2 . 

3. Extrapolate the primitive variables to the intermediate time 
level [a: ] n+ i by homogeneous extrapolation. 

4. Calculate the corresponding state [F°] { 2 and spatial 

source term X i 2 . 

5 . Extrapolate these unknowns from the cell centre [x ' ]; to the 
interfaces at [a; ]^_ i and^ 1 ]^! by stationary extrapolation. 

The resulting unknowns are denoted by u ._ \ andw. /_ 
respectively. 

6. Calculate the states [F°] n _l and [F°] n x 2 _ immediately 
right and left of the cell's interfaces. 

7. At each interface (e.g. at [x } i _i) solve the (approximate) 
Riemann problem defined by the states immediately left 
and right of that interface using the Roe solver: given the 

primitive variables [u] . \ and [u] . 2 and the metric at 

n+i 

the interface [g] . \, determine the characteristic speeds \ k 

l ~ 2 

and the eigenvectors e k at each interface. Project the state 
difference on the eigenvectors: 



8. Update the intermediate state according to 

Ax° 



(8.1) 



[F°] 



o n ™+^ 



" L J * 2 Ax l 

^ [{1 - (?k +((?k - Vk)4>k} hakek]™ + i 

+ ^ [{1 + c k - (cr k - v k )ip k } \ k a k e k ]™*\ 

+Ax°(xr 2+ -x: 



(8.2) 



This completes the recipe for the state update. For every newly 
determined state or flux, new unknowns should be determined. 
After stationary /homogeneous extrapolation the new density is 
available and consequently determining the new unknowns is 
trivial. However after changes in the state other than those due 
to extrapolation, the unknowns must be determined by a more 
expensive procedure (see Sect. 5). The total amount of work per 
cell and per time step in one spatial direction consists of one 
time extrapolation, two space extrapolations, one approximate 
Riemann problem, and one primitive variables determination 
from the state. 

The implementation is such that, by virtue of the station- 
ary extrapolation from the cell centre to the interface, a one- 
dimensional stationary flow in a time independent metric is 



preserved. As the Roe solver only yields the interface flux, and 
the primitive variables that need to be known for stationary ex- 
trapolation cannot always be determined uniquely from a flux, 
extrapolation from the interface to the cell centre is avoided. 
This is the main justification for the seemingly unbalanced pro- 
cedure of applying stationary extrapolation from the cell centre 
to the cell interface but not vice versa. See Appendix C for a 
more elaborate discussion. 

Although the source difference term in step (8) is formally 
required to make the method second order accurate in time, 
numerical experience shows that in most cases its omission 
does not cause a significantly different result. Nevertheless it 
was retained in the calculation of the results presented in this 
paper. 

In case there is no source term and no dependence on known 
functions the same formalism as in Sect. 4 is automatically re- 
covered, because the extrapolations become trivial. Only steps 
(1), (7), and (8) of the recipe are then necessary. 

In Appendix C several alternative implementations are dis- 
cussed and the reasons for choosing this particular one are 
explained. 

For smooth flow the method is second order accurate both 
in space and in time. This is obvious if one realizes that the 
recipe above consists of a second order accurate integration in 
time and a second order accurate integration in space combined 
in such a way as to yield a second order result, see Sect. 3. 
Alternatively one may prove that the methods are second order 
accurate by writing out the nested dependencies. This is done 
in Appendix D. 

9. Flux Limiter 

In this section the algorithm that switches between constant and 
linear interpolation of the cell averaged data is discussed. It was 
developed by Roe for the classical case, and we apply it in GRH 
as well. Because of the different merits of constant and linear 
interpolation it is desirable to switch between Eqs. (4.9) and 
(4.10) depending on what seems the best interpretation of the 
fluid state. Near shocks the jump of conserved variables should 
be reproduced. Although the true solution may be discontinuous 
at any point in space, the limited information available per 
cell generally forces one to place the jump at a cell boundary. 
Therefore near shocks the average state can be considered to 
be a characteristic value throughout a cell, whereas in smooth 
regions of the flow it is better to linearly interpolate between 
adjacent cells. The general interpolation method can be written 
as(seeEq.[4.14]) 



[F i r 2 



u \F'ty:+iF i ty 



2Z^ 



[{u k -(a k -iy k )ip k }\ k a k e k ]^ + £ . (9.1) 



The actual choice between interpolation assumptions is made 
by a "flux limiter" function ip k = ip( r k) of a parameter r k that 
should indicate how smooth or discontinuous the flow is. No- 
tice the assumption, implicit in the notation, that for switching 
purposes the different characteristics can be treated separately. 
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A convenient way to define the parameter r^ is 

r i r 



rk 



n+j- 



[a-kX 



(9.2) 



This parameter is approximately equal to unity in regions of 
smooth flow; consequently, tp(l) = 1 should hold as may be 
deduced from Eq. (4.1 1). Numerical experience has shown that 
"Superbee" is a good choice (Roe 1985): 



( r < 0; 

tp(r) = I min(l,2r) < r < 1; 
I min(2, r) 1 < r. 

In RH the slightly modified form 

r r < 0; 

tp(r) = I min(l,2r) < r < 1; 
I min(1.75, r) 1 < r, 



(9.3) 



(9.4) 



is more robust and yields results that are similar to those ob- 
tained with Superbee (see Sect. 12). Both functions are shown 
in Fig. 4. 



Fig. 4. Flux limiter "Superbee". The original and modified "Superbee" 
interpolation functions are shown as a function of r, the ratio of the 
difference at upwind and the current interface. If the flux limiter is 1 
the data is linearly interpolated between to cells; if it is the upwind 
data is chosen 

This form of the flux limiter is not easy to use numerically 
because the switching parameter is determined by a division 
(possibly by zero). Of course the final update is well defined, 

n+i- 

because of the multiplication by [aj;]. / of the flux limiter. A 
form of Eq. (9.4) that automatically chooses the correct limit is 
found by replacing the combination btp(a/b) with the equivalent 

btp(a/b) =max(0, min(2a, max(6, min(a, 1.756)))) 

+ min(0, max(2a, min(6, max(a, 1.756)))). (9.5) 

Another possible definition of rk is 



rk 



ih] 



[XkO-k]- 



(9.6) 



Indeed the different definitions of r k generally yield comparable 
results except near sonic points. There A& changes sign while 
a,k does not. Thus the second definition of rk yields tpk = 0, i.e. 
shock-like interpolation. This results in unphysical expansion 
shocks, which allow the entropy to decrease along a streamline. 
It is well known that expansion shocks are a problem in first 
order methods of the type described above. Several suggestions 
for a remedy have appeared in literature (e.g. Roe 1985; Van 
Leer et al. 1989; Roe 1992). If the first definition is used, lin- 
ear interpolation is automatically chosen near sonic points and 
smooth transsonic flow is possible and does in fact occur in 
numerical experiments. Thus the first definition of rk is to be 
preferred. 

10. Roe solver for the equations of general relativistic fluid 
dynamics 

The method outlined in Sect. 4 can be used to develop a Roe 
solver for GRH. The vector of all entries of the metric takes the 
place of the known functions g. The vector F° takes the place 
of the state vector W. 

To derive the appropriate matrix A it is convenient to use a 
parameter vector w as in the CNR case (Roe 1981b): 

w = (w a , w ) = (w , w ,w ,w , w ), 



ph ' 



(10.1) 
(10.2) 



where 

w a = Ku a , 

and 

K 2 = ^^ph = -g aP w a w' 3 . (10.3) 

The flow variables can then be expressed in terms of this w by 

r 



K 



„4 



1 



w w ,w w p + Kw g 



4ap 



(10.4) 



By repeatedly using Eq. (4. 16) one can write for two neigh- 
bouring cells 



Fg-F£ = AF a = A a Aw = A a (w R - w L ). 
The Roe matrix of the first spatial coordinate is 

A = A 1 (A )- 1 . 



(10.5) 



(10.6) 



The matrices of the other coordinates follow from symmetry. 
Notice that Eq. (10.6) implicitly assumes that both parame- 
ter vectors are at the same point in spacetime, because terms 
containing metric differences are not allowed for. 
Using the new symbols 

(K)' 
v a = gapv 13 , 



1 + 



(10.7) 
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the A-matrices can be written as follows: 



(K) 



and 





/ C_ — V°Vq 


— V°Vl 






2v V - V *gM Vo 


-v 4 g°°vi 






v 1 - v 4 g m v 


v°-v 4 g 0l vi 






v 2 - v 4 g° 2 v 


-A°V 




\ 3 4 03 

\ v — v g vq 


-v 4 g m v x 




-v°v 2 


-v°V3 -r^T^° 


\ 


-v 4 g°°v 2 


v 4 g m v 3 g 0Q 




-v 4 g m v 2 


v 4 g 01 v 3 «/» 




o _ v 4 g *2 V2 _ 


v 4 g 02 v 3 g 02 






-v 4 g 03 v 2 v° 


- v 4 g Q3 v 3 g m 


) 



(10.8) 



two determinations, in general v a v a ^ — 1 ! The corresponding 
eigenvectors are 



ei 



&2 ■ 



C_ , V 



,1a-. 



„0a„,l 



-(g v - g v ) 

y 



.la„.0 



„0a„,l 



c-,v + -{g v - g v ) 

y 



e 3 = c_ + 



e 4 = 



es = 



r-i' 

c+v 2 ,0,0, 1,0 
c + ^3, 0,0,0, 1 



(10.15) 



A 1 = IK) 



— VVq c_ — V V\ 

v l - v 4 g 01 v v° - v 4 g 0l v l 



4 11 
4 12 
4 13 



2v l — v 4 g n v\ 

2 4 12 

v — v g vi 

3 4 03 

v — v g vi 



—v l v 2 
-v 4 g m v 2 

4 11 
1 4 12 

v — v g v 2 

4 13 



—v l v 3 
-v 4 g Ql v 3 

4 11 
4 12 
1 4 13 

^ — v 9 v 3 



o 13 / 



(10.9) 



and similar expressions for A 2 and A 3 . 

For determining the right eigenvalues A& and eigenvectors 
ej; of the matrix A one can use the relation 



The decomposition of a vector d = (S, A , A 1 , A 2 , A 3 ) in 
terms of these eigenvectors is given by d = J2 o-k&k, where 



ai = 
a 2 = 
a 3 = 



2es 2 
-1 

2es 2 



s 2 k + sy^A 1 - v 1 A ) + (r - l)e(5 + c + v a A a )) 



s 2 k - sy(v°A l - v l A°) + (T- l)e(6 + c + v a A a )) 
(2s 2 k + (T- l)e(6 + c + v a A a )) 



a A = A 2 +- ((g 02 v l - g l2 v°)(v°A l - v l A°) - kv 2 ) 
e 

a 5 = A' + - (G/V - g 13 v°)(v°A l - v l A°) - kv 3 ) , (10.16) 
k = g°°v l A l - g m (v°A l + v l A ) + g u v°A () . (10.17) 



(A l -\ k A°)(A°)-'e k 



0. 



The eigenvalues are then found to be 

Ai=A", 

A 2 = A + , 
A3 = A4 = A5 = A , 



The numerical ease of this decomposition is related to the 
(10.10) determinant of the set of eigenvectors. This determinant is 
— 2es 3 /[(r — l)y]. Hence decomposition is numerically most 
difficult for low Mach numbers. 

To decompose the state difference one should take 



(10.11) 



: = [ F °] T : 



[f°Z2 



(10.18) 



where 



2„01 



± _ (1 - IVKV - s l g m ± sy 



(1 - rVyV> - s 2 g 



„2„00 



A -- 



with 



; Tv\l-v a v a )-UT-l)(l + v a v a ), 



(10.12) 
(10.13) 



This completes the GR Roe solver. The formulae that are needed 
for application rather than derivation are: (10.15) , (10.16) and 
(10.18) to find the decomposition of the state difference (or 
equivalently the flux difference) along the characteristics. 



11. Limiting cases 



and 

,,2 



Before discussing test cases, it is useful and straightforward to 

deduce the Roe solver for a number of important limiting cases. 

It shown in Appendix E that in the limit of smooth flow the 
(10. 14) characteristic velocities found by the Roe solver are indeed the 

characteristics of one-dimensional GRH. Below the Roe solver 
The quantity s is the sound speed in the comoving inertial for a number of special situations is derived by considering the 
frame. Notice that since all these quantities are averages of appropriate limit of the GR Roe solver. 



r E(i-rA + W-A"), 

e = g a W-2g i) W+g n v i) v l) . 
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11.1. Special relativistic Roe solver 

To derive a Roe solver for the equations of special relativistic 
fluid dynamics, all one needs to do is to replace the general 



a/3 



metric g 
9 aP =9c 

r^ 7 = o, 



by the Minkowski metric rfP : 



(11.1) 



This 


yields 






^0 = 


-*°, 






V J = 


v J , 






e = 


v%» - 


V 


V 



k = v u A u -v l A 



i Al 



(11.2) 



Apart from the substitution of the metric, the expressions for the 
eigenvectors and the coefficients cannot be simplified further. 
The expressions thus obtained are similar to the ones discussed 
in Marti et al. 1991 and Marquina et al. 1992. 

11.2. Non relativistic Roe solver in arbitrary spatial 
coordinates 

In order to obtain the non relativistic limit, consider a weak 
gravitational field $ << 1, in a Newtonian coordinate frame 
(X a ). The metric is (Misner et al. 1973, 18.4): 



-(1 + 2$)dX°dX + (1 - 2$)6 ij dX i dX 1 , 



(11.3) 



where $ is the gravitational potential normalized to zero at 
infinity. 

The gravitational potential in the spatial part of the metric 
does not influence the motion of non relativistic fluids because 
it and its contribution to the Christoffel symbols enter only in 
terms that may be neglected in the non relativistic limit. Thus, 
for a frame with arbitrary spatial coordinates and Newtonian 
time (x° = X°, x 1 ) one has 

, „^c dX k dX l 
9oo = -1 -2$6 k r 



dx° dx ' 



goi 



9ij 



dX k dX l 
dx % dx° ' 
dX k dX l 
dx % dxi 



(11.4) 



The spatial part of the four-metric is the three-metric that 
leaves distance invariant. If this is a frame that moves non 
relativistically with respect to a Newtonian frame, i.e. goi << 
^fgTi, one obtains — g — > det(gij) = g c . 

One arrives at the non relativistic limit if in these coordinates 
both the mean and the random velocities of the particles are 
small compared to the speed of light: 

P 

P 

gijU l u J << 1 



«1, 



As a result, 

K 2 ^ Jg~ c p = K 2 . 

Define the non relativistic enthalpy by 

h c = \ \giju % v? + (1 + goo)] + gmu 1 + ■ 



IP 



(11.6) 



(11.7) 



Then 



\ + h c 



K c v? 



V 



T 



i0i 



w A - 


-> K c - = w c , 
P 


rpOO _ 


-+ P, 


rpiO 


-+ pu\ 


rpij _ 


-+ pu l u J +pg v 



1 + -! 

(K e \ 



K c 



\+h r 



(w 



4\ 



K c 



(11.8) 



The speed of sound is 

s 2 -H- (r - 1) {h c - [goivl + \{3ijv\v{ + 1 + ffoo)] } = s\. 

(11 9) 

Notice that both the enthalpy and the sound speed are three- 
scalars that depend only on the cell choice through averaging, 
but not on the three-metric itself. With this sound speed the 
eigenvalues become 



A 






yi Sc 



(11.10) 



Substitution of these limits in the eigenvectors yields, to first 

order, a set of linearly dependent eigenvectors. This reflects the 

well-known fact that the total energy is dominated by the rest 

(11.5) mass energy in the classical limit. Because the latter cannot 
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/frOO 



iin 



classically be converted into other forms, it can be subtracted field by choosing $ = and/or a fluid with no inertial accel- 

from the total energy: erations by taking goi = and/or look at the equations in a 

Cartesian frame. These limits are easily obtained by inserting 

pu — > pn c — p, (1111) the appropriate metric. The expressions for the case of spherical 

pu 1 — > ph c + pg \ coordinates and no gravity are used in Mellema et al. (1991). 

_ . , . , . . , The simplest case, i.e., the combination of all three simplifica- 

This amounts to subtracting the first component of the eigen- . ™ Tr , T T r .l.v ..l ^ ^ 

, ° , r . . tions, is CNRH. In that limit the expressions are exactly those 

vectors from the second. The conserved vector quantity is now „ ,. „ , inoi ,. ... ,. , . „ . XT . 

n J found by Roe (1981b), which were discussed in Sect. 4. Notice 

(11 12) that this is more than just a verification that the Cartesian limit 
of the formulae is correct: the correspondence also hinges on the 
choice of parameter vectors. These were chosen in such a way 
that the non relativistic Cartesian limit is a linear combination 

*• ' ' of the elements of the parameter vector used by Roe: 



T m ' = phu u u u +pg 
f 0i =phu°u i +pg (H 



F° = VTc (p, ph c ~ P, pu 1 ) ■ 

The flux vector in the first direction is 

01 



fg~c (pu , ph c u +pg , pu u 1 + pg 



while the NR limit of the source vector is given by 

s = ^ (o, s a - P (r« + 2r« «•) - r?.r>) , 

and the NR limit of the eigenvectors is 



(11.14) 



ei 



e-2 



e3 



e 4 



es 



l.ftc 



l,h,+ 



=K+s ul ),< 



(v l c + g 01 ),vl + 



1> 2^ 9kivk ^ + l + 900 ^ + 9okV c ' ^ 



0,^20,0,1,0 

0,v 3c , 0,0,1 



Let A = A — S denote the difference in the energy equa- 
tion after subtraction of the rest mass energy. The coefficients 
of the projection reduce to 

r- 1 



a\ 



2s? 



{ \ [g tJ vy c - (1 + fl0 o)] 6 + A - v je & } 



2s cV ^ TT 

r- l 



a 2 



2s\ 

A 1 -v^S 
+ 2 ScV ^ TT 

r- l 



{ \ [ gij vivi - (1 + goo )] 6 + A°- v jc Ai } 



a 3 



a\ 



\(K 



90iV c 

.12 



g ij v l c v> c )& + v jc A : > -A 



A oi 



A z 



v 2 5 - ^(A 1 - v\8), 



,11 



K c [\+h c 



IP 



- v 1 v 2 v 3 £ 

) v ci c ' c ' 



(11.17) 



77.3. Motion in a spherically symmetric gravitational 
potential 

To describe spherically symmetric motion of a fluid under the 
influence of a spherically symmetric gravitational potential, it is 
convenient to choose spherical coordinates (a; = t, x l = r,x 2 = 
(11.15) ^i x3 = < t > ') tnat are related to Cartesian coordinates (X a ) by 



X°=t, 

X 1 = r sin 6 cos < 
X 2 = r sin 8 sinq 
X 3 =rcos6. 



(11.18) 



The only non-zero elements of the covariant metric tensor are 

goo= -l -2$, 
011 = 1, 



gi2 = r , 
333 = r 2 sin 2 ( 



(11.19) 



a 13 , and the volume element is ^fg~ c = r 2 sin#. Note that in fact 

a 5 ^ A - v c S - — (A -v c 6). (11.16) the non relativistic limit has fifn = 1 - 2$ but the gravitational 

potential and the extra non-zero Christoffel symbols it induces 
Three further simplifications of these formulae are possi- are too small to make a contribution to the NR limit. The non- 
ble. One can consider a fluid with no background gravitational zero Christoffel symbols are 



<9$ 



J-oo 




dt ' 


iO r 
01 _ l 10 
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<9$ 
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r 1 

x 22 
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-r, 


r 1 

1 33 
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—r sin 2 6, 


a r 2 

12 _ v 21 
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The characteristic speeds are 

Ai = v\ - s c , 
A 2 = v\ + s c , 
h = v l c , 

corresponding to the eigenvectors 
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(11.25) 



r? 3 = — sin 9 cos 8, 



r 3 -r 3 

1 13 - 1 31 



->3 _ t3 



1 



r^ 3 = r^ = cotg 



(11.20) 



A spherically symmetric solution to these equations is found 
if the source vector and the initial- and boundary conditions are 
spherically symmetric (u 2 = w 3 =0). Then the only non-zero 
components of the stress-energy tensor are 

f 00 = ph c -p, 

f m =ph c u\ 

T = pu , 

T u = pu l u l + p, 

22 P 



T 33 = P_ 



(11.21) 
r~ sur u 

This yields, after dividing out sin 6, the non-trivial vectors 

F° = r 2 (p, ph c - p, pu 1 ) , 

F l =r 2 (pu l ,ph c u l ,pu l u l +p) , (11.22) 

and the source term according to Eq. (1 1.14) 

S = r 0,S V - p— -2pu —,S +2- - p— .(11.23) 
\ at or r or J 

The quantities for numerical approximation become 
K 2 = r 2 p, 

h c = -$ + iwV + ^— --, 

w l = K c u\ 



4 it P 
w =K C -, 



s 2 c = (T-l){h c -\vlvl+$>) 



e\ = (l,h c - s c v c ,v c - s c ) , 
e 2 = {\,h c + s c v\,v\ + s c ) , 
e 3 ={l,±vy e -*,vl). 

The projection coefficients are then given by 



(11.26) 



r- l 



A 1 -^ 



«i = -^~r {[\v]v] + $]5 + A - v\ A 1 ) ^ 



«2 = -^~r (\\v\v\ + $]5 + A - v\ A 1 ) + 



2s c 



r- l 



a 3 = —y- {[k - v\v\]5 + v\ A 1 - A 



(11.27) 



(11.24) 



12. Test problems 

In this section the performance of the general relativistic Roe 
solver is tested on problems to which an exact solution is known. 
These test problems are listed in HSW1 and were used by these 
authors to test five numerical schemes (HSW2). For the precise 
description of all test problems and their analytical solutions we 
refer to these references. The results of our GR Roe solver can 
thus be compared to both the analytical solution of HSW1 and 
the numerical results of HSW2. HSW usually quote the relative 
errors after averaging over the relevant region. These errors 
are not sensitive to zone to zone oscillations. Therefore we 
present the maximum relative errors. With either error definition 
regions where the analytical value of a quantity is zero have to 
be excluded. A measure of the absolute errors there may be 
obtained from the figures. 

Notice that both the test problems and the analytical so- 
lutions themselves are dissipation free. The artificial viscosity 
that explicitly enters some numerical schemes of HSW is only 
a numerical tool to damp unwanted oscillations near shocks. 
Although all numerical schemes have a certain degree of nu- 
merical viscosity, no explicitly introduced viscosity is needed 
to obtain almost oscillation free shocks in the present scheme. 
In all test problems there is no external source term (S a = 0) 
and no time dependence of the metric. The same scheme based 
on the GR Roe solver was used on all test problems. Hence, the 
GR algebraic extrapolation routines, for instance, were called 
even in test problems that could be solved within the framework 
of CNRH. Such problems thus provide an essential test of the 
NR limit of the GR scheme. 
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Table 1. The maximal relative errors of the Roe solver and several HSW schemes for the Riemann shock tube (test problem 1). 
HS W errors have been estimated from their figures 



method 


At 


Ax 1 


i u ° 


u 1 


e 


V 


nhti. u " 


"V-s 00 


ptlUi ^/- g m 


Roe 


25 


1 


0.75% 


0.21% 


0.75% 


0.27% 


0.76% 


HSW Wilson 


25 


1 


7% 


8% 


7% 


7% 


14% 


HSW Barton 


25 


1 


2% 


2% 


3% 


2% 


2% 


HSW Mono 


25 


1 


4% 


3% 


2% 


4% 


2% 


Roe 


12.5 


0.5 


0.58% 


0.18% 


0.58% 


0.24% 


0.58% 


HSW Wilson 


12.5 


0.5 


7% 


8% 


7% 


7% 


11% 


HSW Barton 


12.5 


0.5 


2% 


2% 


3% 


2% 


1% 


HSW Mono 


12.5 


0.5 


3% 


3% 


2% 


3% 


2% 



Table 2. The maximal relative errors of the Roe solver and several HSW schemes for the Riemann strong shock tube (test problem 
2). HSW errors have been estimated from their figures 



method 


At 


Ax 1 


•, u " 




e 


V 


phut ," 


r V~s m 


Roe 

HSW Mono 


0.75 
0.75 


1 
1 


6.62% 
12% 


1.45% 
5% 


8.60% 

21% 


3.19% 

25% 


6.01% 
17% 


Roe 

HSW Mono 


0.15 
0.15 


0.2 

0.2 


7.24% 
6% 


0.28% 

2% 


6.72% 
10% 


1.12% 
12% 


1.69% 
11% 



Fig. 5. Comparison of the numerical (crosses) and analytical solution of the Riemann shock problem on a 100 zone grid (test problem 1) 



Fig. 6. Comparison of the numerical (crosses) and analytical solution of the Riemann shock problem on a 200 zone grid (test problem 1) 



Fig. 7. Comparison of the numerical (crosses) and analytical solution of the Riemann strong shock problem on a 100 zone grid (test problem 2) 



Fig. 8. Comparison of the numerical (crosses) and analytical solution of the Riemann strong shock problem on a 500 zone grid (test problem 2) 

12.1. Non relativistic Riemann shock tube tests 

Test problem 1 is a non relativistic shock in a gas with T = 7/5. 
The physical domain stretches between x l = and x l = 100. 
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Across a membrane at x l = 35, which is removed at a; = 0, the 
gas is characterized by 



(pu Q ) L = 1 10 5 , 

4 = 0, 

PL = 1, 



(p U % = 1.25 10 4 , 
u l R = Q, 
PR = 0.1. 



(12.1) 



The numerical solution is calculated on a 100 zone grid 
with Ax 1 = 1 and a time step of A a; = 25. The solution at 
x° = 5000 is compared to the analytical solution in Fig. 5. 

The contact discontinuity is spread over 4 cells instead of 
the 6 cells in their best solution and at the same time the zone 
to zone oscillations are far less pronounced. The solution is 
already so accurate that little improvement may be expected 
from further grid refinement. For the sake of comparison the 
above test problem is also solved on a 200 zone Aai ' =0.5 grid 
with half the time step used above: Ax° = 12.5. The solution, 
again at a; = 5000, is plotted in Fig. 6. Compared to the solution 
on the 100 zone grid there are indeed hardly any differences. 
The errors are shown in Table 1 . 

The Roe solver produces solutions that are almost an order 
of magnitude more accurate than the best solutions of HSW. 
The solutions of all schemes improve relatively little on a finer 
grid. This is an indication that the errors are dominated by the 
discontinuities, which are scale free. 

Test problem 2 considers a strong non relativistic shock 
problem in a gas with T = 5/3. It is specified by: 



(pu°) L = 1 10 2 , (pu°) R = l, 



,i - 



= 0, 

PL =2/3, 



,i - 



= 0, 

PR = 2/3 10- 1 . 



(12.2) 



The numerical solutionis calculated on a 100 zone grid with 
Aa; 1 = 1 and a time step of Aa; = 0.75. In this problem the 
method gives a negative pressure in cell 38 after 4 time steps if 
the flux limiter as suggested by Roe is used. A negative pres- 
sure is lethal to a characteristic-based scheme since the speed 
of sound cannot be determined. By using the slightly modified 
form of Sect. 9 no problems arise. If one applies the original 
Superbee on the coefficients of the flux difference (A^a^) also 
no negative pressures arise, but with this modification expan- 
sion shocks form instead of transsonic expansion waves. The 
solution for a; = 172.5 is compared to the analytical solution in 
Fig. 7. Results on a 500 zone grid with Aa; = 0.15, are shown 
in Fig. 8. 

The performance on this test problem shows the ability of 
the relativistic Roe scheme to accurately form and propagate 
NR shocks, contact discontinuities and rarefaction waves. The 
errors are given in Table 2. Again the Roe solver compares 
favourably to the HSW scheme. This time the contact discon- 
tinuity and the shock are hardly separated on the coarse grid 
and consequently the accuracy is improved on the finer grid. 
Nevertheless the results of the Roe solver on a 100 zone grid 
are comparable to (if not better than) the HSW solution on a 
500 zone grid. 



12.2. Special relativistic Riemann shock tube tests 

Next, four special relativistic problems were treated. A rela- 
tivistically moving gas with T = 5/3 collides at the centre of 
the grid (a; 1 = 50) at a; = with a similar gas moving at an 
equal speed in the opposite direction. This problem is com- 
pletely equivalent to the wall shock problem of HSW, except 
that it is presented in a boundary independent way. The initial 
condition is 



(pu") L = l, 



(PU°)R = 1, 



u l L = \/vPvP- 
V l = 2/3 10" 



u R = -Vw°w°- 1, 
p fl = 2/3 10- 6 . 



(12.3) 



First, we take u° = 2.24. The numerical scheme was run on 
a 100 zone Aa; 1 = 1 grid with a time step of Aa; = 0.5. The 
result is shown in Fig. 9. The most prominent error is the spike 
in internal energy at the point of collision. This "wall heating" 
phenomenon is the actual solution of the finite difference equa- 
tions as was shown by Norman & Winkler (1986), and appears 
in all methods tested by HSW as well. Apart from this spike, the 
result is excellent compared to the analytical solution. Contrary 
to the solutions of HSW the post shock density and specific 
internal energy do not have a significant systematic error. In 
addition the zone to zone oscillations are far less pronounced 
than with the HSW scheme. The errors are given in Table 3. 

The Roe errors are typically a factor three less than those 
of HSW. The methods of HSW become unacceptably unstable 
around u° = 4. They reach the same conclusion as Norman 
& Winkler (1986, presented in 1983), namely that extremely 
relativistic shocks could -at that time- only be reproduced by 
implicit codes, which are much more expensive than explicit 
methods. The present explicit method, however, performs well 
on extremely relativistic shocks: results for u° = 625 and T = 
5/3 are shown in Fig. 10 and for u° = 625 and T = 4/3 in 
Fig. 11. 

The errors are given in Table 4 and are in fact somewhat 
smaller than those at u° = 2.24. This is caused by the fact 
that the Courant number of these calculations is closer to unity, 
which has a favourable effect. 

Recently, Marquina et al. (1992) have solved this test prob- 
lem up to u° = 70 with a somewhat different implementation 
of a Roe solver. Hence, the Roe method is to our knowledge 
the first explicit scheme that can deal with extremely relativis- 
tic shocks. This is strong evidence that the advantages of the 
characteristic based schemes for violent NR flows (Woodward 
& Colella 1984) are retained by the relativistic version of the 
Roe solver. 

In test problem 4 a special relativistic shock in a gas with 
T = 5/3 is considered. It is specified by a membrane at a; 1 =45 
and: 
(A = 10, (pu°) R = l, 



<=0, 
PL =40/3, 



u R = 0, 

p R = 2/3 10~ 7 . 



(12.4) 



The numerical solution is calculated on a 100 zone grid 
wifhAa; 1 = 1 and a time step of Aa; = 0.5. The solutionis for 
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Table 3. The maximal relative errors of the Roe solver and several HSW schemes for the mildly relativistic (u° = 2.24) collision 
(test problem 3a). HSW errors have been estimated from their figures 



method 


At 


Ax 1 


•, u ° 


u° 


e 


V 


phut ," 


r V-s m 


Roe 

HSW Mono 


0.5 
0.5 


1 
1 


8.72% 
20% 


0.00% 
0% 


9.58% 

22% 


0.15% 

3% 


0.00% 

4% 



Table 4. The maximal relative errors of the Roe solver and several HSW schemes for the ultra relativistic (u° = 625) collision 
(test problems 3b and 3c). HSW schemes, applied to this problem, are unacceptably unstable 



method 


At 


Ax 1 


•, u ° 


u 1 


e 


V 


phut ," 


r V-s°° 


Roe (r = 5/3) 
Roe (r = 4/3) 


0.5 
0.5 


1 
1 


4.16% 
4.99% 


0.00% 
0.00% 


4.34% 
5.25% 


0.01% 
0.21% 


0.00% 
0.00% 



Fig. 9. Comparison of the numerical (crosses) and analytical solution of the mildly relativistic (u° = 2.24) collision (test problem 3a) on a 100 
zone grid 



Fig. 10. Comparison of the numerical (crosses) and analytical solution of the ultra relativistic (u° = 625) collision with T = 5/3 (test problem 
3b) on a 100 zone grid 



Fig. 11. Comparison of the numerical (crosses) and analytical solution of the ultra relativistic (u° = 625) collision with T = 4/3 (test problem 
3c) on a 100 zone grid 



Fig. 12. Comparison of the numerical (crosses) and analytical solution of the relativistic shock tube (test problem 4) on a 100 zone grid 



Fig. 13. Comparison of the numerical (crosses) and analytical solution of the relativistic shock tube (test problem 4) on a 500 zone grid 

x° = 50 is compared to the analytical solution in Fig. 12. The errors in Table 5 show that the Roe solver results on the 100 
solution on a 500 zone grid with Ax 1 = 0.2 and a time step of zone grid are better than the HSW results on a 500 zone grid. 
Ax° = 0.1 is compared to the analytical solution in Fig. 13. The The improvement in the Roe results on the finer grid is again 
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Table 5. The maximal relative errors of the Roe solver and several HSW schemes for the relativistic shock tube (test problem 4). 
No HSW results are available on a 100 zone grid. HSW errors have been estimated from their figures 



method 


At 


Ax 1 


i u " 


u 1 


e 


V 


nhti. u ° 


r V-s m 


ptlUi ^/- g m 


Roe 


0.5 


1 


3.00% 


0.69% 


1.51% 


2.71% 


1.94% 


Roe 

HSW Mono 


0.1 
0.1 


0.2 

0.2 


0.04% 

7% 


0.66% 

3% 


1.05% 
14% 


2.58% 
8% 


0.11% 
13% 



due to the fact that the contact discontinuity and the shock were 
not well separated on the coarse mesh. 

12.3. General relativistic tests 

In test problems 5 and 6 the performance of the code in a 
general relativistic setting is tested. Both problem are situated 
in a Schwarzschild metric. The region of a; 1 = 2 to x l = 18 is 
covered by 16 equidistant cells. 

In problem 5 radial geometric infall of dust is treated. The 
analytical solution is again discussed by HSW. Unfortunately 
their last expression for the analytic value of the energy (HSW1 : 
Eq. [67]) is incorrect. The correct expressions are, in the present 
notation (a; 1 = r): 

d(x 1 )= -£hL = - 



(x 1 )- 



e{x l ) = 



Vix 1 ): 



1\ _ _E_ 



r-i 



K,p r /(r-i) 



(z') : 



(12.5) 



'1 



As an example the test shown by HSW, with boundary 
condition T = 5/3, D = 1 .6 10~ 2 huo = — 1 and k = 12/23 was 
run. The initial condition of the gas is given by D = 1.6 10~ 4 
huo = — 1 and k = 12/23. The problem was solved with two 
different time steps Ax° = 0.1875 and A a; = 0.5. The result 
at x° = 90 is shown in Fig. 14. The errors are summarized in 
Table 6. 

The errors in the dynamically dominant density and velocity 
are extremely small (< 0.1%). The relative error in the energy 
is still small, but large compared to the other relative errors 
(3%). This is due to the fact that in the conservative form of the 
equations the total energy is solved for. Hence the error in the 
internal energy is dominated by small relative small errors in 
the much larger density and velocity. In the (non conservative) 
formulation of HSW the rest mass energy is subtracted from the 
total energy and then the internal energy is solved for. This is 
numerically more benign in case of a dynamically insignificant 
pressure. The errors on the coarser grid are less than on the 
finer grid. This is explained by the fact that the fluid is to evolve 
to a steady state, which analytically is only reached after an 
infinite amount of time, due to time dilatation near the black hole 
horizon. Numerical diffusion, which is larger on the coarse grid, 



accelerates this evolution. This is not the explanation for the 
HSW errors, because they have reached a numerical steady state 
at the time the calculation is stopped. If given sufficient time 
the errors of the present scheme become arbitrarily small by 
virtue of the algebraic stationary extrapolation subgrid model. 

In this problem the use of algebraic extrapolation really 
pays off: implementations based on numerical integration of 
the source term (see Appendix C) succeed in reproducing the 
dynamically dominant density and velocity but fail to repro- 
duce the (dynamically insignificant) pressure, due to numerical 
diffusion. 

To test the ability of the present method to deal with 
transsonic accretion onto a black hole, a Bondi-Hoyle accretion 
problem was solved (test problem 6). This differs from geomet- 
ric infall in the sense that the pressure is no longer negligible. 
The pressure support results in a lower accretion velocity. The 
analytical solution is found by assuming a radius for the sonic 
point x s , a value for the radial component of the four- velocity 
w' and for the entropy constant k. Given the critical radius, 
the metric at the sonic point is known: g 00 = —1/(1 — 2/x l s ), 
g 11 = 1 — 2/x\. Upon substitution of the metric and the as- 
sumed u\, Eq. (6.23) yields h/ho. The latter constant and the 
requirement of sonic flow combine in Eq. (6.25) to yield the 
enthalpy h. Using the assumed k, the sonic point density p s 
may be found from Eq. (6.20). The pressure at the sonic point is 
found from Eq. (6.10). With all these variables known the two 
remaining constants of the flow (D and huo) are easily found. 
The Bondi-Hoyle problem HSW show is given by: T = 5/3, 
x\ = 8, u\ = 1/4 and k = 120/23. This yields h s = 26/23, 



Ps = 
huo 



1 10 3 and flow constants D 

= h s (u ) s 



{<YKp 



1.6 10- 



■26/23 -y/13/16. It may be verified that these 
values are indeed a solution of Eq. (6.1 1). As an initial condi- 
tion the analytic solution for D = 1.6 10~ 4 , huo = — 1 and 
k = 12/23 was used. 

The problem was run to a; = 360. Two different time steps 
were employed Aa: = 0.1875 (identical to HSW) and Aa: = 
0.5 (to demonstrate the method's performance with a larger 
Courant number). The results are shown in Fig. 15 and Table 7. 
The errors are at least an order of magnitude less than those 
of HSW. This should again mainly be attributed to the use of 
stationary extrapolation. The main test in this problem is finding 
the sonic point. 

The two-dimensional flows that HSW show were not pur- 
sued here because they are especially geared towards accretion 
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Table 6. The maximal relative errors of the Roe solver and several HSW schemes for general relativistic dust accretion (test 
problem 5). HSW errors are as mentioned in their text. HSW refine the grid towards the horizon; the average spacing has been 
given here 



Roe 

HSW Wilson 

HSW Mono 



0.1875 


1 


0.1875 


0.67 


0.1875 


0.67 



0.09% 
3% 

1% 



0.04% 


2.87% 


1% 


30% 


2% 


10% 



2.93% 
9 



method 


At 


Ax 1 


•, u ° 


u 1 


ir u ° 


V 


phut ," 


r V-s m 


r "V-s m 


Roe 


0.5 


1 


0.06% 


0.01% 


0.68% 


0.69% 


0.04% 



0.06% 
9 



Table 7. The maximal relative errors of the Roe solver and several HSW schemes for general relativistic Bondi-Hoyle accretion 
(test problem 6). HSW errors are as mentioned in their text. HSW refine the grid towards the horizon; the average spacing has 
been given here 



method 


At 


As; 1 


•, «° 


u° 


n C U " 


V 
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Fig. 14. Comparison of the numerical (crosses) and analytical solution for general relativistic dust accretion (test problem 5) on a 16 zone grid 



Fig. 15. Comparison of the numerical (crosses) and analytical solution for general relativistic Bondi-Hoyle accretion (test problem 6) on a 16 
zone grid 



flows, hardly allow quantitative comparisons and yet require 
substantial amounts of computational effort. 



13. Discussion 

In this paper we have presented a method for the numerical ap- 
proximation of solutions of the equations of general relativistic 
fluid dynamics. The method is based on Roe's approximate Rie- 
mann solver for the classical Euler equations. A disadvantage 
of the method is that it is rather cumbersome and thus computa- 
tionally expensive. However, it is physically appealing because 
of the direct physical meaning of characteristics as the paths 
along which information propagates. The economics is thus 
determined by whether the extra effort results in a sufficiently 
more accurate scheme. If so, the flow can sufficiently well be 



approximated on a coarser grid. This coarser grid pays back the 
extra work per cell, which becomes especially advantageous as 
calculations are done in more spatial dimensions. As shocks are 
quite difficult to reproduce numerically, it should be expected 
that the more violent flows benefit the most from the above 
scheme. Indeed, there is evidence from the above test problems 
that shows that the Roe scheme can be applied on a significantly 
coarser grid than the HSW schemes to obtain similarly accurate 
results. In addition the extremely relativistic shocks such as the 
ones that occur in the test problem of colliding special relativis- 
tic fluids (see Section 12), that until now required an implicit 
treatment, can be handled accurately by the explicit Roe solver. 

The method presented in this paper has been successfully 
applied in two-dimensional flow simulations of interacting 
winds. It was discovered that a wide range of parameters pro- 



Frits Eulderink, Garrelt Mellema: General relativistic hydrodynamics with a Roe solver 



25 



duces hydrodynamically collimated jets. The collimation mech- 
anism was called "inertial confinement", because it is a mech- 
anism by which jets form when a spherical hot wind blows 
outward through a thick accretion disk. Results for non rela- 
tivistic planetary nebulae are given by Mellema et al. (1991), 
Icke et al. (1992), Mellema (1993), Mellema (1994), Mellema 
& Frank (1995); relativistic jets were calculated by Eulderink 
& Mellema (1994). 

The code has not been yet applied in the context of active 
nuclei. But in cases of violent discontinuous flows we expect 
this to be worthwhile given the performance on the test prob- 
lems. 

As in Cartesian non relativistic hydrodynamics (Roe & 
Baines 1982; Roe 1986b), operator splitting to obtain quasi 
one-dimensional equations is an area of possible improvement 
for shock capturing schemes. In view of the above results it 
seems worthwhile to try to extend work on genuinely multi- 
dimensional characteristic-based methods, once proven to be 
successful in Cartesian non relativistic hydrodynamics, to the 
relativistic domain as well. The same holds for alternative im- 
plementations of the Roe solver which eliminate the possibility 
of unphysical fluxes. 



14. Conclusions 

1. The translation of Roe's concepts for an approximate Rie- 
mann solver from the non relativistic to the general rela- 
tivistic domain is laborious but straightforward. 

2. The implementation in a Minkowski metric completely 
analogous to the implementation in the Cartesian non rela- 
tivistic hydrodynamics case presented by Roe (1986a), can 
be accomplished and is successful. 

3. The implementation in GR or for curvilinear coordinates 
can be accomplished via a non-trivial subgrid model to deal 
with the source terms. 

4. A source term splitting exists that allows algebraic integra- 
tion of the curvature source terms. 

5. The numerical scheme built around the GR Roe solver 

yields excellent approximations to one-dimensional NR 
and SR shock tube problems, 

is an explicit numerical scheme that overcomes the prob- 
lems experienced with some other explicit schemes, 
when dealing with extremely relativistic shocks (ex- 
amples up to u° = 625 have been given), 
produces excellent approximations to the steady solu- 
tion of GR radial accretion problems. 
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Appendix A: Determination of the primitive variables from 
the state 



In this appendix three methods to determine the primitive 
variables given the conserved quantity vector are described: an 
analytical method, a one-dimensional Newton-Raphson method 
and a six-dimensional Newton-Kantorovich method. 

A..1. Analytical method 

The first two ways of solving for the primitive variables amount 
to finding the root of a quartic polynomial. That polynomial 
results after substituting the mass and energy components of 
the state in terms of the primitive variables in the relation that 
is found if velocity normalization is used to replace the three 
momentum components by one relation: 
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Then C > 1 because u Q u° > -g 00 , and < C < 1 be- 
cause pressure is positive. Therefore the coefficient of the fourth 
power is non-negative. It may also be deduced that physically 
allowed values for £ are in the range 



/2-r 

0< \l— — <£< 1. 



(A5) 



Proof that precisely two real roots for all allowed coefficients 
exist, and that only one of these roots obeys the physical con- 
straint that it is positive, is presented first. 
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The roots of the polynomial Eq. ( A. 1 ) can also be expressed where 
as the intersection points of two polynomials: 



(ci - \)e i 




70 = 4a l(a 2 - a 4 S) + [(a 2 - a\8) 2 - Aa^a^ 

71 = \ U + ^r? 2 - 4«) , 



i\ -1 



5 = 



1 



V\V 2 



/?30l+^2) + 



ii 

a 1 + <jt + t 2 



2/3|(<7 + T + \a 2 ) 



(A 10) 



1 + 



-C 




-C 



The auxiliary symbols used are 
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The right-hand side polynomial is negative everywhere ex- 
cept on the interval between its roots. Thus if Co = 1 the 
existence of precisely one intersection in the physically accept- 
able region is easily proven. If Co > 1, the polynomial on the 
left-hand side has one root that is smaller than zero, a double 
root at £ = and one root £ > 0. This polynomial is positive 
everywhere except for the region between its two roots. From 
Eq. (A. 6) it may now easily be deduced that there is one root in 
the physically acceptable region 
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(A. I) Q nce £ j s known, one finds by combining the expressions 
for C, Co, T 00 and using hC£, = 1 that the time component of 
the four-velocity is 
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The density is now easy to determine: 



Thus precisely one physical solution exists of the quartic 
Eq.(A.l). 

In the classical limit C j 1 and Co | 1, which causes 
«4 to vanish. The two roots that remain are £ — > 1 and £ — > 
— 1 /(2r — 1). The physical constraints imply that the first root is 
the ratio sought, because the latter is smaller than zero. The root 
formula presented below does indeed pick this correct solution 
in the classical limit. Because the root of the polynomial is 
given below as a continuous function of its coefficients and 
the coefficients in their turn are continuous functions of the 
conserved variables, this means that the proper root is selected 
in the relativistic case also since such a root was proven to exist. 

Having selected the one root that satisfies the physical con- 
straints it can be determined by the exact algebraic quartic root 
formula. However, special care must be taken, because the co- 
efficients of the fourth and third powers of the polynomial tend 
to zero in the NR limit. The familiar expressions for the root of 
a quartic equation should thus be rewritten in a form that can 
accommodate this limit behaviour. 

The selected root £ is found to be 
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The enthalpy, pressure and velocities are then determined by 
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A..2. One-dimensional Newton-Raphson method 

The above procedure is computationally rather expensive be- 
cause of the many extra variables and exponentiations involved, 
and because special care that must be taken when numerically 
evaluating the cubic roots. An alternative way to find the ratio 
is based on an iterative procedure. The obvious choice is to use 
Newton's method on the polynomial Eq. (A.4) . 

With some effort it may be proven that the derivative of the 
polynomial of Eq. (A.4) is always at least 2/T in its physical 
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root, thus guaranteeing quadratic convergence of the Newton 
method. This is the main reason to choose £ and not some other 
variable. With the starting value that is asymptotically correct 
in the NR limit (i.e. as C, C — > 1) 



6 



= 1, 



(A 15) 



only five iterations of the Newton method where required to 
obtain a relative accuracy of 10~ 16 in all test problems. This is 
a lot less expensive than the algebraic formula presented above. 

A.. 3. Six-dimensional Newton-Kantorovich method 

The third alternative to solve for the primitive variables is to 
use a Newton-Kantorovich method on a system of equations. 
This method is presented below in a form that yields 
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Simple explicit expressions exist for the primitive variables 
in terms of K and the elements of the parameter vector. To 
determine K and the parameter vector given the state vector of 
a relativistic fluid one must solve a coupled set of equations. 
The first equation gives a relation between the variables, and 
the other five determine the relation of the parameter vector to 
the state: 
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where the determinant is 
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These increases determine a new parameter vector. Because 
Eq. (A. 19) is in fact non-linear these steps must be iterated to 
obtain convergence. 

In principle one can substitute the first relation of Eq. ( A. 19) 
in the other five to eliminate K. The reason this is not done 
and K is treated as an independent variable instead is that 
(A 19) one circumvents the evaluation of the square root. This is not 
only faster but also avoids problems that otherwise arise if 
convergence goes through a region of negative K 2 . 



Given a starting value for the parameter vector, solve for the 
increases in the elements of the parameter vector that would 
make the right-hand side equal to if the equations where 
linear in the chosen variables. These increases are given by 
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Appendix B: Example of a super-resolution problem 

In this appendix we will give a specific example of how a false 
assumption about the subgrid model can influence the flow. 
Consider NR stationary radial outflow with supersonic velocity 
m 1 . The NR limit of stationary extrapolation in the 1 -direction 
yields that the mass flux 
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the transverse velocity 
u k , k = 2,3 

and the entropy 
V 



(5.3) 



(5.4) 



are invariant. In the NR limit the density p is to be solved from 
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In spherical coordinates (radius r, polar angle 8 and az- 
imuthal angle <j>), the density at radius r is thus determined 

by 
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Now consider the same flow at the same point in space in 
cylindrical coordinates (indicated by ~): r = r sin 8, z = r cos 8 
and 4> = <j>. The total velocity is split in a polar component 
u l = m 1 sin 8 and a component parallel to the axis of symmetry 
u 2 = u l cos 8. Extrapolation to a different r is done in two steps. 
A step in the 2-direction that leaves all quantities the same and 
a step in the 1 -direction that yields 
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Apart from the different pre-factor there is an even more se- 
rious problem with these formulae: there is a pole for sonic flow. 
In the spherical coordinates where the velocity is supersonic the 
density decreases away from the centre. In cylindrical coordi- 
nates however, the subgrid model has a supersonic or subsonic 
1 -velocity depending on the projection of the total velocity, 
which in its turn depends on the position on the grid via 8. In 
case it is subsonic the density increases outward from the axis. 
This is not just a sign error (which is bad enough). Because of the 



pole the difference is extreme at the position on the grid where 
the projected velocity is near sonic, i.e. sin 8 = l/M. However, 
as the flux varies smoothly through M = 1 it may not be affected 
quite so badly. Nevertheless, in principle instabilities can result 
from the false assumption of the subgrid model that stationary 
spherical flow can be approximated by stationary cylindrical 
flow. In numerical experiments on multi-dimensional flows one 
should be aware of these possible instabilities. 



Appendix C: Implementations for general relativistic flow 

In this appendix first the method that results if operator splitting 
as discussed in Sect. 3 is combined with an arbitrary extrapola- 
tion of Sect. 6, is described. Next, it is shown how to arrive at the 
method presented in Sect. 8. Functions that are known a-priori 
are collectively called "known functions". The "unknowns" are 
functions of a minimal set that needs to be known in addition to 
the known functions to completely describe the state and that 
need to be determined by integration of the partial differential 
equations. In the GRH setting one may think of the primitive 
variables. 



C..1. Basic method 

Straight application of operator splitting as discussed in Sect. 3 
and any source term splitting as discussed in Sections 3 and 6 
results in the following implementation: 

1. Given the state [F ]? determine the unknowns u™, 

2. Extrapolate the unknowns to the intermediate time level 
[x°] n+ i_ assuming 



F° Q = T. 



(CI) 



Call the unknowns u i 2 . For this extrapolation the un- 
knowns u change with time as 



u fi = (FX l (T-F"g fi ) 



(C2) 



3. Extrapolate the unknowns at cell centre [a: 1 ]; to the inter- 
faces at [a; 1 ]^! and [a; 1 ]^!, assuming 



F\ = X. 



(C.3) 



The resulting unknowns are called u. \ and u. ? re- 
spectively. For this extrapolation the unknowns u change 
with position as 



u x ={Fl)-\X-Flg A ). 



(CA) 



4. At each interface (e.g. at [a; 1 ]^!) solve the (approximate) 
Riemann problem defined by the local known functions and 
the unknowns immediately left and right of that interface. 

This yields the unknowns [u] . \ and [u] . \ . 
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5. Using Eq. (C.4) , extrapolate the unknowns [it]. \ and 

% 2 + 

n+— + n+— + 

[it], 2 from the interface to the cell centre to form [wL_ 2 

»+ 2 

and [u] i+ 2 , respectively. 

6. Determine the fluxes that belong to these unknowns: 

i n+i+ i n+i+ 

[F%_ 2 and[F'] i+ 2 . 

7. Update the state according to 



[f»]T~ 1+ = [fV 



^{IF^ + -IF^ + 
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8. Given the state [F°] i 2 determine the unknowns u i 2 . 

9. Using Eq. (C.2) , extrapolate these unknowns in time, to 
form the final update [F°]f +1 . 

Since the update sequence begins and ends with an integra- 
tion over time, these integrations may be combined to yield 

[.F ]; 2 immediately from [F°] { 2 . Thus one can dispense 
with steps (1) and (9) if one is only after a next update. 

In case there is no source term and no dependence on known 
functions the old formulae are automatically recovered, because 
the extrapolations become trivial and only steps (4), (7), and (8) 
are then necessary. As is shown in Appendix D, the extrapo- 
lation in time must be at least second order accurate, while 
the extrapolation in space only has to be first order accurate. 
This difference in required accuracy arises because in the spa- 
tial direction information can be obtained from neighbours on 
both sides whereas in the time direction one has only the past 
to predict the future. More specifically, the second order spa- 
tial extrapolation terms cancel in the flux difference. In general 
extrapolation serves to determine the unknowns at a new loca- 
tion in spacetime. Extrapolations that yield combinations of the 
known and unknown functions must therefore be followed by 
a step in which the unknowns are determined. In any case, a 
procedure as described in Appendix A is necessary to determine 
the unknowns after the changes to the state of step (7). 

The total amount of work per cell and per time step in one 
spatial direction is: 

1 time extrapolation 

4 space extrapolations 

5 unknowns determinations after extrapolation 

1 (approximate) Riemann problem 

2 unknowns determinations after the Riemann solver 
1 unknowns determination from the state. 

For smooth flow this implementation is second order accu- 
rate both in space and in time. This is obvious if one realizes 
that the recipe above consists of a second order accurate inte- 
gration in time and a second order accurate integration in space 
combined in such a way as to yield a second order result. Al- 
ternatively one may prove that the methods are second order 
accurate by writing out the nested dependencies. This is done 
in Appendix D. 



C..2. Algebraic extrapolation 

A special case of the above scheme results if the algebraic sub- 
grid model of Sect. 6 is chosen. After stationary /homogeneous 
extrapolation the local density is available and consequently 
determining the new unknowns is trivial. The use of analytical 
solutions for extrapolation makes the GRH scheme as robust 
as the NRC scheme. It allows the use of the same methods that 
guarantee a physical prediction in CH to guarantee an physi- 
cal prediction in the GRH case. More specifically, if one uses 
algebraic extrapolation, a flux limiter that produces physical 
interface states for the Riemann solver and a Riemann solver 
that yields a physical interface state (such as the exact solver), a 
physical approximation for the updated state is guaranteed even 
in GRH! This property is due to two facts: first, the extrapola- 
tion and the Riemann solver predict physically allowed states. 
Second, only vectors that exist at the same place and time are 
combined to yield an update. Notice that if the extrapolation 
equations are integrated numerically, numerical inaccuracies 
can still cause the result to be unphysical. 

Advantages of this implementation are: one-dimensional 
stationary solutions are automatically recognized as such, spe- 
cific use is made of the coordinate symmetries, and the concept 
of operator splitting in the coordinate directions is strictly ap- 
plied. Hence, changes in the state due to changes in the source 
term after the passage of a discontinuity are consistently taken 
into account. The operator splitting results in equations to which 
physical solutions exist. 

The total work involved is 

1 time extrapolation 
4 space extrapolations 

unknowns determinations after extrapolation 

1 (approximate) Riemann problem 

2 unknowns determinations after the Riemann solver 
1 unknowns determination from the state. 



C..3. Determination of the unknown variables from the flux 

The use of extrapolation that requires knowledge of the un- 
knowns such as in the previous subsection in combination with 
a Roe solver has a disadvantage: a Roe solver only yields an in- 
terface flux and not the interface state or the interface unknowns 
(see Sect. 4). Stationary extrapolation of the flux without knowl- 
edge of the primitive variables fails because only the constants 
D and Dhu e can easily be determined from a given flux, but not 
the entropy. Yet, this entropy is required to find the cell centre 
flux. 

The unknowns can be determined uniquely from the state, 
but in general not from the flux. If 



D 2 + ga6(V=gT la )(^T l ?) < o, 



(C.6) 



there exist two different entropies that are consistent with a 
given flux. This is related to the possibility of a stationary 
shock: the entropy may change while the flux remains constant. 
Thus supersonic and subsonic primitive variables may exist that 
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yield the same flux. Different branches may have be chosen on 
either side of an interface, because of the possibility of a steady 
shock. If the mass flux D = the flow is even less determined: 
w 1 = but the three transverse constants hu e and the entropy 
are only related through the pressure at the place where the flux 
is known. 

There are two ways to resolve this problem. One can either 
use additional information to choose a branch or one can try to 
circumvent having to evaluate the unknowns from the flux. In 
this subsection the former path is taken and in the remaining 
subsections of this appendix the latter path is explored. 

In view of the non-uniqueness one must use extra informa- 
tion, such as the initial states that entered the Riemann problem, 
to determine unknowns from fluxes. The unknowns can be dif- 
ferent on either side of the interface due to the presence of a 
stationary contact discontinuity or shock. Therefore a physical 
solution to the mathematical problem above is to use the en- 
tropy and transverse constants determined by the last available 
state on either side of the cell interface. 

A possible way to determine an interface state instead of an 
interface flux by a (first order) Roe type method is 
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where Yl an d Yl denote summation over all indices of 
negative and positive characteristic speeds, respectively. 

Different interface states are predicted immediately left and 
right of the interface, because there exists no single interface 
state if the states immediately left and right are separated by a 
shock or contact discontinuity that has no velocity with respect 
to the grid. Rather there are two distinct states immediately left 
and right of an interface in such a situation. These states cannot 
be used directly to predict the interface flux, because the two 
approximate interface states do not necessarily correspond to 
the same flux. Hence, such a scheme would not be conservative. 
Yet, they can be used to find on each side the branch of the 
primitive variables that correspond to the Roe flux. Notice that 
waves with characteristic speed zero do not affect the interface 
states. Picking the wrong entropy branch near sonic points does 
not cause large errors in the flux. 

C..4. Time extrapolation only 

Several alternatives exist to prevent extrapolation of the inter- 
face flux. To minimize the amount of work one may be tempted 
to choose T = S and X = 0. This corresponds to a subgrid 
model that prescribes a constant flux. Spatial extrapolation is 
thus trivial. After the extrapolation of cell interface fluxes to 
the cell centre no unknowns need to be determined because the 
fluxes suffice. Hence the total amount of work is 

1 time extrapolation 

space extrapolations 

3 unknowns determinations after extrapolation 

1 (approximate) Riemann problem 



unknowns determinations after the Riemann solver 

1 unknowns determination from the state. 

However, there are also drawbacks. One-dimensional sta- 
tionary solutions are not automatically recognized as such. Un- 
physical interface quantities are likely to occur if the metric is 
space dependent. No specific use of the coordinate symmetries 
is made. After the addition of the flux difference an interface 
state must be evaluated while curvature terms have not yet been 
evaluated, which again could lead to unphysical intermediate 
results, especially if discontinuities cross the cell. 

In this case the unknowns that enter the Riemann problem 
must be determined after the extrapolation of the flux from 
the cell centre to the cell interface. In this case the uniqueness 
issues may be solved by selecting a branch after numerically 
integrating the extrapolation equation for the unknowns. 



C..5. Centre-centre extrapolation 

Another way to avoid the determination of the primitive vari- 
ables from the flux for solvers that only yield a flux is to com- 
bine the extrapolation before and after the approximate Rie- 
mann problem to just an extrapolation before the approximate 
Riemann problem: the unknowns at the centre of cell i — 1 
are extrapolated to the centre of cell i. There these unknowns 
and the unknowns that correspond to the centre of cell i are 
combined by an (approximate) Riemann solver to yield a flux. 
Likewise a flux can be formed by the Riemann solver based on 
the extrapolation of the unknowns of cell centre i + 1 and the 
unknowns at [a: 1 ];. The two fluxes determine the update of the 
state in cell i. Because the above extrapolation is not symmetric, 
two Riemann problems must be solved per interface instead of 
one. 

At first one might be tempted to consider the approximate 
Riemann problems as Riemann problems at [x ];, because that 
is where the unknowns exist. However to make the scheme 
second order accurate the transport speeds that enter the recipe 
via the Courant number must be evaluated at the interface as 
must the eigenvectors. Also, a scheme that takes these quantities 
at the cell centre is not conservative for equations without a 
source term: the eigenvalues and/or the eigenvectors of systems 
of equations depend on the state (and possibly the metric). Thus 
the projection of the same flux on the local eigenvectors is in 
general different from one place to another. Consequently, a 
part of the flux that is considered upwind at [a: 1 ]; may not be 
considered downwind at [a: 1 ];-!. Therefore the flux that enters 
cell i may not be the same as the flux leaving cell i — 1 . To resolve 
this problem one should project the above flux difference on the 
eigenvectors at the interface, although the fluxes are defined 
at the centre of cell i and have no physical meaning on the 
interface. To obtain the unknowns at the interface the cell centre 
unknowns must be extrapolated to the interface as well as to the 
cell centre of their neighbour's cell centre. Using the interface 
characteristics guarantees that if there is no source term and 
consequently the flux difference vector at [a: 1 ];- 1 is the same as 
the flux difference vector at [a: 1 ];, the approximations to their 
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state differences come out the same too. Thus the calculated 
interface flux is the same in those cases. This guarantees a 
conservative scheme if no source terms are present. 

However, if there is a source term in one or more but not 
in all equations of a set, the two flux differences above are 
no longer the same (as a vector), and it is possible that the 
conserved quantities corresponding to the equations without 
source term are not conserved by this scheme. This may not 
be as bad as it sounds, since the proper treatment of shocks is 
the strongest motive for placing emphasis on conservation and 
shocks are likely to be treated well by the method; nevertheless, 
a conservative scheme is preferable 

C..6. Approximate extrapolation 

To save the extra work involved in the algebraic extrapolations 
or to take external source terms into account, one may consider 
numerical integration of the source term. Numerical integration 
of the spatial extrapolation equation is straightforward as it 
only needs to be integrated to first order accuracy. The time 
extrapolation has to be second order. The scheme below is 
especially apt to deal with numerical extrapolation: 

1. Given the state [F°] { 2 determine the unknowns u { 2 . 

2. Calculate the source term X i 2 . 

3. Extrapolate the unknowns to the time level [x°] n 

M? = Mr' + + l 2 Ax° [{F°)-\T - F,Vo)]"- i+ • (C8) 

4. Extrapolate the unknowns to the intermediate time level 
[x°] n+ i via 

[ut +h =Mr' + + A a; [(^)- 1 (T-^Vo)];. (C.9) 

5. Calculate the corresponding state [F°] i 2 and spatial 

source term X i . 

6. Extrapolate the unknowns from the cell centre at [a: 1 ]; to 
the interfaces at [a: 1 ]^! and [a; 1 ]^!, assuming 

[<|; = Mr 1- - ^ mr\x - Fi g , x) ]f 2 - , 

(CIO) 

7. Calculate the states [F°] n _l and [F°] n x 2 _ immediately 
right and left of the interfaces of each cell. 

8. At each interface (e.g. at [a; } i _i) solve the (approximate) 
Riemann problem defined by the local known functions and 
the unknowns immediately left and right of that interface. 

This yields the flux [F 1 ]^ . 
% 2 

9. Extrapolate the fluxes to the cell centre, based on the old 

source term: 

[F l ]Z k2+ = [F i ri+ l 2 Ax l [X]T 2 -, 



10. Update the state via 

[FQ] n + i + = [F0] n + i 



Ax° 



+ Aa; ([X]""' + -[X] n r 2 ~ 



[F% 



(C12) 



ift+- 2+ =if%i-\^ i ixt: 



(Cll) 



For this method the extrapolations are so cheap that they 
are ignored. The amount of work is: 

time extrapolations 

space extrapolations 

3 unknowns determinations after extrapolation 

1 (approximate) Riemann problem 

unknowns determinations after the Riemann solver 

1 unknowns determination from the state. 

This method involves less work than the alternatives men- 
tioned above, but it has the disadvantage that one-dimensional 
stationary solutions are no longer recognized and that the lin- 
ear extrapolation in space and time may not yield physically 
allowed states. Particularly a discontinuity may cause the inter- 
mediate result to be unphysical. 

In the update only the spatial flux differences enter; the 
fluxes themselves do not. Therefore for a second order scheme 
the estimates for the fluxes immediately left and right of an in- 
terface are sufficiently well approximated by first order extrap- 
olation from the cell centre, since the second order terms cancel 
in the flux difference. This allows one to use first order extrapo- 
lation. An implementation based on approximate extrapolation 
was used to obtain the results of Mellema et al. (1991), Icke 
et al. (1992), Mellema (1993), Eulderink & Mellema (1994), 
Mellema (1994), and Mellema & Frank (1995). 

C..7. Mixed extrapolation 

The method that was used in this paper is based on the source 
term splitting of Sect. 6 and on algebraic extrapolation from the 
cell centre to the interface and linear extrapolation followed by a 
source term correction in the opposite direction. This may seem 
unbalanced but it offers the advantage that a one-dimensional 
stationary solution is recognized, while at the same time the 
difficult step of determining unknowns from the Roe flux is 
evaded. The amount of work is then 

1 time extrapolations 

2 space extrapolations 

unknowns determinations after extrapolation 

1 (approximate) Riemann problem 

unknowns determinations after the Riemann solver 

1 unknowns determination from the state. 

Time extrapolation is trivial in all examples mentioned in 
this paper, because none of them have a time dependent metric. 

Appendix D: Proof of second order accuracy 

In this appendix it is proven that for smooth flow the method 
presented in Appendix C is second order accurate both in space 
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and time, by writing out all the nested dependencies. Only Likewise one has 

terms that are relevant to prove that the method is second order . 

accurate are kept. ' [X] n *f w \ ( [Xf + _X~_ + [Xf + _\ 



The update for the state is A o 

[F°q + = [ft"- - §j£ <yq + - t^q + ) - « i (^f-1" + r*r 1_ ) 

+ IA,°[X.(^)- 1 (X-^)]^- 

[F°q- « [F ]? + lA^mJ 1 ~ 2 { [X]? -' + 2 Ak0 fe.O +*«(*£)" ' (^ - ^Vo)]".! 

+ A(Az ) 2 [T 3ff ,o + TJ^)" 1 (T - F° fl ,o)] , + M + 5 A *° \ X >9fi + ^(^V (^ - ^V<C} 

+ IA,°[X.(^)- 1 (X-^)]^-, 

L ^q + - ^Tt^A^txr!*, [X q + « i f^q- + [X ti 



2 

{ Ax° 



2 Ax 

if% + > ■ « [j^q- ' - iA^txr- ' , w i / [Jf r j- + [X]r; 



^(^Virr [^rr - [f 1 ]" 



il n n + 2+ _ rZ7h n + 2 + 1a lry 1 n+ 2 + 



[^'tj =i|[ fl ]3- + [- F Vi^ _ 1 J IVI „ , 1 Ajfv „ „,.„,,-l/ f71 „,! 



„0 



+ IA,°[X.(^)- 1 (X-^)]^ 



2 



> 1 ]^-+[F 1 ]T ' 



+ IA,°[X U (^)- 1 (X-F 1 1 )]; + ?"- (^-2) 

1 2 

+ 2 Ax [F U (F U ) [X — F l )\ i _\ After substitution of these expressions one finds 

i { [F 1 ]?., + lA^ [F>, + i^V (T - F,Vo)]r_, [FT'" * t^r'" + A," [(X - Fl)]" + |(Az ) 2 { 



[T]" +5+ « [T]" +I + A/i;^ )-' (* - F\) 



[FTr^UiFTr+iF^ 



[TXl + \Ax a [T g g iQ + T u (F°y l (T - F°g fi )] 



1 Ac ! ! n+i- / ! n+i- rE ,l n ™ + 



A* 1 u u «+i V * + i+ *,-7 + As°T u (F°)- 1 (*-.F , 1 1 i). (I>.3) 

~ i /'rc l h n+ 2 — r7?h n+ 2 — ^ 

~ 2 ^L-r J; + L-c* Ji+i J After substitution of these expressions in Eq. (D.l) one can 

, conclude that the terms of first and second order in A a; are 

+ 2 Ak [F u (F u )~ (J-Fi)] , 2 indeed equal to their analytical counterparts 

Ax°[S - F\yy (DA) 

and 

+ [F^ +i + IA,° [F>, + i^V (T - ^, )]; +1 } i (Ak o )2 { ^ q + ^ (i ,o r i (5 _ j,! _ ^ 



[F 1 ]? + ^A* [F>, + i^ )" 1 (T - F° fli0 )]" +1 
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respectively. 

Appendix E: Characteristic form of the equations 

In this appendix it is shown that the present method applied 
in the limit of two (almost) identical interface states makes 
use of the exact characteristic form of the equations of general 
relativistic fluid dynamics projected onto one spatial dimension. 
For two identical states the averaging in Eq. ( 1 0.7) becomes 
trivial: 

< K > ^(y^gph)\ 



P_ 
ph ' 

ph 



(E.l) 



Substitution in the other expressions of Sect. 1 is now straight- 
forward. 

The hardest to check is the expression for the sound speed 
in arbitrary coordinates x a , given that in Lorentz-Minkowski 
coordinates X a which are comoving with the instantaneous 
systematic fluid velocity the sound speed is s. Because the fluid 
is at rest in this frame, 



dx a dx? dx c 



dr dX? dr dX ' 
Thus in terms of the systematic velocity u° 



dX i dXi 



(E.l) 



(E3) 



Having established these relations, maximize the sound speed 
in the 1 -direction over the spatial components of the rest frame 
sound speed, 



A=^ 



ax 1 

81° 



5° 



8x" On 

ax" D 



(EA) 



under the conditions that the rest frame sound speed is s: 



S° = 



1 



VT 



VlJ s*si = - — - 2 . 

1 — s z 



(E.5) 



Using a Lagrange multiplier p one obtains the equations for the 
spatial components of the rest frame sound speed 



dx° 
A^r— = pS : , j = 1,2,3. 



(E.6) 



dXi dXi 

These relations can be substituted into the restrictions 
Eq. (E.5) to yield, after replacing the products of partial deriva- 
tives by the metric elements: 

2 * 2 



1 



= (g 00 + u°u°)X 2 - (g 01 + uW)2X + (g u + u l u l ). 



(E.l) 



With this relation for the Lagrange multiplier applied to the 
spatial components of the rest frame sound speed, one obtains 

(1 - s 2 )(u°u°X 2 - 2u°u l X + u l u l ) = s 2 (g m X 2 - 2g m A + g 11 ). 

(E.S) 

The roots of this equation are exactly equal to the roots X^ of 
Eq. (10.1 1) in the limit of small differences. 
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